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A GENERAL ALGORITHM FOR THE 
CONSTRUCTION OP CONTOUR PLOTS 


Wayne Johnson 

Ames Research Center 
and 

Aeromechanics Laboratory^ . 4.^-4 
AVRADCOM Research and Technology Laboratories 


Fred Silva 
Informatics, Inc. 


SUMMARY 


An nKjoritta U de«:rlbed that nerfcm, the task o£ draeinq equal- 
level contovre on a plane, which requites interpolation in two 
dimensions based on data prescribed at points distributed irreg- 
ularly over the plane, •rhc approach is described in detail. The 
computer program that implements the algorithm is documented and 

listed. 
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1.0 INTRODUCTION 

The qraphical presentation of experimentally or 
theoretically qenorated data sets frequently involves the con- 
struction of contour plots. Consider a dependent variable z that 
is a function of two independent variables x and y; z = f(x,y). 
The functional form f is not known. It is assumed that f is a 
sinqle-valuod function of x and y. By measurements or calcula- 
tions, the value of z is obtained at a set of N discrete points. 
The data may be presented in qraphical form in terms of contours 
of equal value of z on the x-y plane. To construct such contours 
it is necessary to interpolate the values of z between the ore- 
scribed data points. In qeneral, these data iK>ints may bo dis- 
tributed irreqularly over the x-y plane. This report describes 
an algorithm developed to construct contour plots for such cases. 
The computer program that implements the algorithm is documented 
and listed. 

1 . 1 Poser ipt ion of the Approach 

The tlatvi arc prescribed at a set of N points distributed 
irregularly over the x-y plane: 2 ^, for n=l to N. In 

order to perform the interpolation, the sxsints on the x-y plane 
are connected by straiviht line segments, to form a set of tri- 
angles with a convex boundary (figure 1). Then the data can bo 




interpolated over the cdqet of the trianqles. To construct the 
contour for the value Z, all points on the edqes where arc 
located. Finally, these points are connected to form the z^7. 
contour, with the trianqulation alqorithm described here, the 
interpolation alomi the edqes often involves widely separated 
points on the x-y plane. In such a case, linear irterpolation 
between the end points of the ed< 7 e is unlikely to nroduce a 
smooth contour. Hence, it is usually necessary to smooth the 
data, by usinq a least-squared-error fit of the data to a bi- 
variable polynominal for z=f(x,y). Then the interpolation 
alonq the edqes is performed usinq this functional form. it is 
also possible to use some standard technioue to draw a smooth 
curve throuqh the interpolated noints on the edqes of the tri- 
anqles. in summary, the alqorithm involves four basic steps: 

(a) trianqulation of the plane? jb) smoothlnq of the data; 

(c) interpolation alonq the edqes; and (d) drawing the contours 


The first step is trianqulation of the plane. There are 
N data points and The tr ianqulation will bo d.oscribod by 

an array that identifies the two end points of each e<lqe, and an 
array that identifies the three vertices of each trianqle. At 
each staqe in the procedure, there are a set of points that 
define (in order) a boundary, inside which the triangles have 
been identified. At the start, all the data points are outside 
the boundary, and no points on the boundary have been located. 
The last data point and the data point closott. to it are used to 



start the procedure: they are the initial boundary points, and 

are no lonqer in the set of points outside the boundary r one 
edge has been identified. Thereafter, the algorithm proceeds 
by marching around the boundary, examining points outside the 
boundary relative to a boundary edge. The objective is to 
identify a point that together with the boundary edge wiil 
form a new triangle. The criteria for locating such a point 
are that it be closest to the boundary edge and that there be 
no other points within the resulting triangle. These criteria 
are satisfied by locating the point such that the oarameter n 
is minimized, where D equals the sum of the distances from the 
point to the two end points of the boundary edge. The points 
examined in this manner are those on the boundary, immediately 
adjacent to the boundary edge being considered; as well as the 
points outside the boundary. From the points outside the 
boundary it is necessary to exclude any for which the resulting 
triangle would overlap the triangles already identified 
(within the boundary) , which requires two tests. First, relative 
to the boundary edge there is a side within the boundary. The 
straight line formed by the boundary edge and its extensions to 
infinity divides the x-y plane into two half-planes. All 
points that are either on this line or in the half-piano corre- 
sponding to within the boundary are immediately excluded. 

Second, the point identified as closest to the boundary edge is 
examined to determine whether the two new edges of the resulting 
triangle would pass through any of the boundary, inside which 
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the trianqles have been identified. At the start, all the data 
points are outside the boundary, and no noints on the boundary 
have been located. The last data noint and the data point 
closest to it are used to start the procedure: they are the 

initial boundary points, and are no lonner in the set of noints 
outside the boundary; one edqc has been identified. Thereafter*, 
the alaorithm proceeds by marchinn around the boundary, oxamininq 
points outside the boundary relative to a boundary edqe. '’’ho 
objective is to identify a noint that toqether with the boundary 
edqe will form a new trianqie. The criteria for locatinq such 
a noint are that it be closest to the boundary edqe and that there 
be no other noints within the rcsultinq trianqie. These criteria 
are satisfied by locatinq the point such that the parameter D is 
minimized, where D equals the sun of the distances from the noint 
to the two end noints of the boundary edqe. The points examined 
in this manner are those on the boundary, immediately adjacent to 
the boundary edqe beinq considered* as well as the points outside 
the boundary. From the noints outside the boundary it is 
necessary to exclude anv for which the resulting trianqie would 
overlap the trianules already identified (within the boundary) , 
which requires two tests. First, relative to the boundary edqe 
there is a side within the boundary. The straiuht line formed 
by the boundary edqe and its extensions to infinity divides the 
x-y Plane into two half-nlanes. All points that are either on 
this line or in the half-plane correspondinq to within the 
boundary are iminediately excluded. Second, the point identified 



as closest to the boundary ed<re Is examined to determine whether 
the two new edges of the resulting triangle would pass through 
any of the other edges on the boundary (which may happen if the 
boundary is concave). If so, the point is excluded. When a 
point has been successfully found from among the points outside 
the boundary, a new triangle and two new edges have been identi- 
fied; a new boundary point is inserted between the two current 
boundary points being considered (hence two new Ijoundary edges 
replace the old edge) ; and the point is no longer outside the 
boundary. When a ooint has been successfully found from among 
the adjacent boundary ooints, a new triangle and one new edge 
has been identified; and the middle boundary point is no longer 
on the boundary (hence the new boundary edge replaced the two 
old edges). This procedure continues, marching around the 
boundary until there are no more ooints outside the boundary. 

The boundary may be concave at this stage, however, so the 
procedure still continues, examining adjacent boundary ooints 
relative to each boundary edge until the boundary is completely 
convex, that completes the triangulation. The end points of all 
edoes have been identified. For the interpolation procedure it 
is necessary then to identify those edges that form the boundary. 
To draw the contours, the four other edges that form the two tri- 
angles on either side of each edge must be identified as well. 
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The followiiwi relationships are useful. Let P «» number 
of data points, F. » number of edqes, T » number ol trianqles, 
and B - number of boundary noints or edqes. Then 

E » jT + jD 

^ * jT + (jB + 1) 

so T » 2(P - 1) - B 

E « 3(P - 1) - B 
e-t = p- i 


The minimum number 
number of trianqles 
The maximum number 


T . = P-2 and E_. 

min min 


of boundary ooints B . = 3 qives the maximum 

^ luin 

and edqes: T^^^^ = 2P-5 and = 3P-6. 

of boundary points is = p, which qives: 

= 2P-3. 


The trianqulation deoends only on the x and y coordinates 
of the data noints, hence it is the same for all dependent var- 
iables. The remaininq steps depend on the deoendent variable as 
well. 

The second steo in the alqorithm is smoothina of the data 
for z. This steo is ootional, and does not depend on the tri- 
anqulation. The 2 -surface is fitted to a noLynomial of the form: 

I L . . 

z = z :: c,. 

i=0 j«0 

' 3 ' 



where 


K » maximum (I>J) 

L ■ minimum (K-i^J) 

The innut oarameters I and J define the highest powers in the 
polynomial. The coefficients c.^ are obtained from a loast- 
souared error fit of this function z to the actual data z, at 
the set of N data points. Then the polynomial is used to 
evaluate a new set of z values at the data point. This set of 
smoothed values of the dependent variable reolaces the original 
data in the interoolation algorithm. The error of the smoothed 
data is defined as: 




N 

n=>l 


(z 


■ z_ 


’old 


) 

’now 


The third step is interpolation alonq the edqes. The 
contour value is specified. Then each edqe is examined to 
determine whether where and z^ are the values of 

the dependent variable at the end points of the edge. If so, 
then there is a point on the edge where z = Z, hence this is a 
Doint on the required contour. This point is obtained by linear 
interix)lation between the end points if the data has not been 
smoothed. If the data has been smoothed, the fitted polynomial 
is used to evaluate z along the edge and hence locate the ooint 
where z = Z. The result of the interpolation procedure is a 
Qf noints on the olane where z ^ Z, and the edges on 



which these points are located. 


The fourth step is drawing the contour for z = 7,. The task 
ts to convert the interpolated points in the proper order. The 
contour will consist of one or more I ines that either start and 
end at a boundarv edge, or are closed curves. There can only 
be one contour through a triangle. The nroccdurr starts by 
searching the list of interpolated points for one that lies on 
a boundary edge. There are two outer edg.-r that form a triangle 
with this boundary edge, which wore idontifie.i in the triangula- 
tion algorithm. The contour must pass through one, and only one 
of these edges. So the list of interpolated points is searched 
for the point that lies on one of these two edges. There are 

four edges (identified in the triangulation algorithm) that form 

two triangles with one edge on which this second point lies. 

The list of interpolated noints is searched for the f)oint that 
lies on one of these four edges. (There will, be only one such 
point in the list: one from each of the two triangles, and one 

of these will be the previous point on the contour.) The nro- 
cedure continues searching for points in this fashion until 
another boundary point is reached. Then a contour lino is 
drawn through these points, in the order located. The procedure 
is repeated until there are no more points in the list that lie 
on boundary edges. If there are atill interpolated points that 
have not been used, there must be a contour segment that forms 
a closed curve. One of the remaining points is nicked as a 
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starting point, and tha abova procadura is followed until this 
starting point is encountarad again. Than a contour line is 
drawn through these points, in the order located. The procedure 
is repeated until all the interoolated points have been used. 

The desired contours are specified in terms of a base value 

z and an increment Az, so the contour value is Z « + nAz 

o 

where n is any integer (positive, negative, or zero). The 
interrelation and contour drawing steps are repeated for every 
such Z that lies within the range of the data. 

The computer program described here does not include the 
graphics software. The user must supply the subroutine that is 
called to draw the contour on the particular graohics device 
being used for the output- 

1.2 Summary of Component Modules 

The above procedures are computationally independent steps 
in the process. For this reason, each procedure is self- 
contained within separate subroutine modules. One master sub- 
routine is called by the user program and it, in turn, controls 
and sequences the execution of the procedures described above. 

The master subroutine also accepts, by means of an argument list, 
the date and parameters that the user supplies for the procedure. 
In addition, the user supplies a subroutine for graphics output 
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the contour Un« « they are generated. 

alqorlth™ for oertarn „,^„ten»ined arid on the 

data points define a redu „ian,uUtlcn sub- 
plane. It may be deeira distribution of 

tisM^Iceuent -iU often increase the execution s.ed 
points. ^ ^ , ,arce nueber 

apbstantiaUy. h o- 

data points given _ triangle edges. 

to allow for a linear „ot be exercised and 

the smooth inq op^ 

For such a case, fi^tinn could be deleted 

d for the surface curve fitting cou 
the procedure substantial s.avinds in cb,ect 

altoqether. This would result 

time program size. 

• i. =.« he used to modify the 
other variations which may 
"" of reducing obiect time storage require 

method for the nurnosc «,dif ications are 

„«„ts or increasinn execution speed 

discussed later in Section 7. 

- of this section is composed of a short 
yne remainder of this 

description of each component module. .. fid 

heirarchy diadram of the orocessind pacKad • 
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CNTLNS 

This is the subroutine accessed by the users callinq proiiram 
Cor drawing contour lines of constant for some set of data 
defining z = f(x,y). CMTI.NS is supplied with the known values 
of x,y and z, several computational parameters, and a list of 
constant z values for which contours arc to be caleulateil and 
drawn. There must be at least I triplets of x-y-z points and 
no duplicate points are allowed. The function z = f(x,y) must 
be single valued. 

TRIANG 

Called by CNTLNS. This subroutine constructs the convex 
polygon of triangles from the x-y data. 


MIDDLE 

function subprogram used by TRIANG. This routine finds the 
middle value of three known integer values. 

S MSRF 

Called by CNTLNS. Performs least-square smoothing of the 
2-surface. The smoothing is an optional procedure. 
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LLSOF 


Called by SMSRF. This Is a utility module taken from the 
International Mathematical and Statistical Library (IMSL) . 

IJiSOF is used to solve a linear least-squares problem. It 
solves for the solution vector X of the general problem 
AX * B, where A is the coefficient matrix and H is the riuht 
hand solution vector. LLSQF is a proprietary program; t.LSOF or 
its equivalent must be obtained by the user. 

INTERP 


Called by CNTLMS. Performs linear or non-linear interpolation 
over the triangle edges for constant contour values. 

PO LYX2 

Function subproaram used by INTFRP to evaluate the polynomials 
obtained in SMSRF for values on trianulo edges. 

CNT OUP 

Called by CNTLMS. Reorders interrelated points into proper 
contour lines, noth closed and open contours are accommodated. 
CNTOIJR calls a user supplied subroutine to draw tli i contour line. 
The user subroutine must bo named CNTCRV, 
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CBVCHK 


Called by CNTLNS. If the user specifies a base value and 
increment scheme for defining 7.^ (as described later) , then 
this routine is used to verify that is within the ranee of 
the known data, if not, Z^ is incremented or decremented by 
AZ until is in the proper ranee. 

CNTCRV 

called by CNTOUR. This is the user suoplied subroutine used to 
draw the contour on the nraphics device. 






2.0 


MASTER SUBROUTINE 


The subroutine CNTLNS is the user's application proqraro 
contact with the contour software, its primary function is to 
check for errors and, based on user incut oarameters, control 
and properly sequence the calls to other modules which perform 
the computational tasks. After all requested contours have been 
processed, control is passed back to the application proqram. 

2.1 Pescription of Argument List 


CALL CNTLNS (X, Y,Z ,N, TSMOPT, IKXP, JEXP.NCNTRS.CLIST, 
EPSLON, TKRR) 


Input arguments; 


Xp = the list of independent variable values for the f motion 
'Z = f(XjV) for n = 1 to N 

'n “ independent variable values for the function 

2 = f(x,y) for n - 1 to N 

”n -■ the list of vlependent variable values for the function 

■^ = f(x,y) for n ^ I to N 

N = the ranqo of N for the x,y and z lists 

ISMOPT •- smoothing option parameter 

= 0 for no smoothing 

0 then the function z = f(x,y) is smoothed by means of 
a least squared error curve fit 

lEXr = highest order of the smoothing polynomial for x if 
ISMOPT 0 
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je:«p 

NCNTRS 


CLISTj 


Return 


EPSLGN 

lERR 


■ highest order of. the smoothing polynomial for Y if 
ISMOPT fi 0 

(The dimension C in the program must be at least 
(K+l) (L+l-JjK) where K ■ min(I«J) and L » max(I,J).) 

■ the number of contours of constant Z to be generated » 
and NCNTRS <50. If NCNTRS < 0, then the program will 
determine constant Z values to process from the relation 

Z * Z nAZ 

CO 

where Z » constant 7. value 

Z^ * contour base value 
o 

AZ * increment value . 

« If 1 < NCNTRS < 50, then CLIST is the list of constant 
Z values (Z ) For which contours will be generated, for 
j»l to ncntSs. 

If NCNTRS < 0, then CLIST(l) is taken to be Z^ and AZ =• 
CLIST (2). “ 


arguments; 


= the error c, introduced by the smoothing if ISMOPT 0. 

» return error flag 
=s 0 for no errors 

= 1 for N<3 or N'MAXPTS where MAXPTS is the maximum number 
of x,y,2 triolets allowed 

- 2 for invalid lEXP and/or JEXP values if ISMOPT 0 
(Note - 

lERR is 2 if the number of coefficients resulting 
from IRXP and JEXP is greater than MAXCOF or greater 
than N, the number of points under consideration) 

(Where MAXPTS is the dimension N, and MAXCOF* is the 
dimension C in the program.) 
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Required dimensions ; 


X(N) 

Y(N) 

Z(N) 

CLTST(50) 

ZNKW(N) 

IE(R,2) 

ITE(E,4) 

XI (E) 

ETA(E) 

LAJ'.BDA(E) 

IDE(E) 

IPOWR(C) 

JPOWR(C) 

COEF(C) 

For the array dimensions given above, and for all array dimensions 

used in this document, the following definitions apply r 

N = the maximum number of data points to be processed 

C • the maximum number of coefficients to be used for 
smoothing 

E = 3N-6 = the maximum number of triangle edges oroduced 
by the triangulation of N points 

T = 2N-5 = the maximum number of triangles oroduced by 
the triangulation of N ooints. 


2,2 Descr ipt ion of Algorithm 


Figures 3a and 3b present a block diagram of the mudule CNTLNS. 
The functions of parts A to M are as follows: 
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Piqufe 3a. Block nianram of CNTLNS, Parts A to P 
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for^errorl variables and check input arguments 

MAXCOP ■ 23, MAXPTS * 500 

lERR - 0 

E » 0.0 

U N<3 or N>MAXPTS qoto 



(B) Call subroutine to triangulate the X-Y data 
CALL TRIANG(X,y,N,LEDCES,IE,IBE,ITE) 

(C) If smoothing option is off (equal to zero) then skip 
around sections D and ^ 

if ISMOPT * 0, goto i^riK 


(D) Else, check the request exponents for errors 

if IEXP<0 or JEXP<0, goto^9^ 
n - iEXP+1, J1 - jExp+iC/ 

NMIN a minimum of II, ji 
NMAX a maximum of II, Ji 

if J1>I1 then MC » (IEXP+1 ) * (jexp+1) - texp/ 5 i 
IT Jl<ll then NC = (JEXP+1)^E.XP+1 - JEXP/2?^ 
iT NON or NOHAXCOF, goto^ JEXP/2) 


{Tor Kal,(^XCGF 
|IP0WR(K) a 0 
|£POWR(K) = 0 


Z-f the function 

Call SMSRF ( X , Y , Z , ZNEW, N , I EXP , JEXP , NCOEP , COEP , I POWR , JPOWR ) 

calculate 

if NCOEF <0 go^/l 2 ^ 


for k a I to N 
(j^/FLOAT (M) 


fllOj for k * 1 to N 

values Which can be a^co^aLS! the contour 
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(ro 

(H) 


ZMIN » 2(1) 

2MAX » ZMIN 

(Tor k * 2 to N 
ZMIN « minimum of ZMIN, Zj^ 
[£max ■ maximum of ZMAX, Zj^ 


Branch around the next section if the contour list is qivon, 
otherwise call subroutine to range check tho base value and 
r' 3 set it if necessary. 



K »0 
FN-1.0 




if NCNTRS * 0 qoto USOy 

call CBVCHK (cSST(WcLIST( 2) ,ZMIN,ZMAX,CLNEW) 
if CLIST(l) ^ CLNEW then CLIST(l) » CLNEW 


(I) 


Determine the (next) contour constant value. 




K 


K+1 


( 1 ) 


ZCON = (K) (FN) (CLIST(2)) + CLISX(1 
if ZCON>ZMIN and--z CON < ZMAX goto pScX 
IT FN^O.O goto HOm 
FN =>-1.0 
K 

goto / 


130) K 



if K>NCNTRS 00^(300) 
ZCON = CLIST(K) V_y 


(J) 

(K) 





ZCON = CLIST(K) 

If ZCON<ZMIN or ZCON^ZMAX goto (20^ 

Call subroutine INTERP to interpolate 
contour points for constant Z 

CALL INTERP ( X , Y , Z NEW , ZCON , LEDGES , IE , lEXP , JEXP , I3MOPT , LAI ’BDA , 
XI , ETA , J , COEF , I rOWR , JPGWR , NCOEF , N ) 



if J;i^0 , CALL CNTOUR { ZCON , XI , ETA .LAMBDA , J , IBE , ITE) 

goto 

RETURN 


lERR = 1 
RETURN 

lERR ® 2 
RETURN 
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2 3 Description of Subrouti ne CBVCHK 

Thl» subroutine is celled by CNTI.S3 after the 7. date ran.ro 
has been determined. CBVCIIK mill check the .liven value of the 
contour base var.-e (V to verify that it is within the ranee 
of the data. If not, the base value is shifted by the diven 
increment (azl until ZHAX, and the shifted vaiue of z„ 

assumes the new reset value. This verification and resettimi 
of % is often useful for cases in which the eiven base value is 
only°a auess by the user and the tance of the 7, data may not be 
known in advance. The areument list (or CBVCIIK is established as 

follows; 


CALL CBVCHK { ZZERO , DELZ , ZMIN , ZMAX , /sZrEW) 


Where and iz are the selected base and increment values for 
selectin, the constant values of Z for the contours, ZUh and 
XMiVk are the data ranee as calculated in CXTLIIS. . is, " 

return, the base value which falls between 7.MIN and ZMAX and may 

or may not be equal to Z^- 
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3.0 


TRIANGULATION SUBROUTINE 


The subroutine TRIAN6 performs the triangulation/ as 
described in Section 1. This subroutine uses the function 
middle. 

3 . 1 Description of Argument List 

CALL TRIANG (XD,YD,N ,L,E,BE,TE) 

The triangulation algorithm is supplied with a set of N 
data points , i=l to N. The coordinate pairs are to be 

connected by straight lines to form the triangles. The 
procedure input consists of: 

XD(i) “ the list of abscissa values 

YD(i) = the list of corresponding ordinates 

N = the range of i; the number of points 

in the x and y lists 

The subroutine output consists of a set of index oointers 
defininq each triangle edge, each boundary edge of the final 
polygon, and indices of adjacent edges to each triangle edge. 
The subroutine output is stored as; 

E(^»2) = index pointer of end points of a triangle 

edge in ascending order (E( i , 1) <E( i ,2) fo 
all «; for i « 1 to L 


2S 



BE(IU 


TE(i,4) 

L 


• 1 if the i-th row of E defines a boundary 
edge; otherwise equal to zero; for 
it St 1 to L • 

« index of adjacent edges for each 

corresponding row of E; for i » i to L. 

« total number of edges constructed by the 
triarigulation procedure 


An assortment of local variables are used during the 
triangulation process and are defined as follows: 


P(j) 


J 

B(k) 

K 

T (m, 3) 


M 

X(i) , Y(i) = 


Index numbers of points lying outside 
the boundary of the triangulated ^ints. 

P lists the indices of the remaining 
candiaate points, for j « 1 to J. 

Number of points remaining in array P. 

Index numbers of points defining the 
current boundary, in order, for k = 1 to K. 

Number of values listed in array B. 

Indices of triangle vertices of each 
triangle, in ascending order, for 
m 1 to M. 

Total number of triangles; the same as 
the limit of m for array T- 

Arrays of X and Y data after the XD and 
YO input values have been scaled by the 
ra. ge of data. Scaling of the data 
eliminates problems with machine precrsion 
while leaving the relative nosi ion of the 
data points unchanged. 


3.2 Desc ription of the Algorithm 

Figures 4a to 4e present a block diagram of the module 
TRIANG. The functions of parts A to Y are as follows. 


26 



initial ize local variables 
and scale the XY data. 


■rake the last X-Y pair 
(i-th pair) as the first 
»x)int . 


rake the noint nearest X j , 
Yi, to be the second point 
of the first ed^ie. 


Store the two rapints as the 
first boundary edqe • 




Fet oointers (Kl and K2) 
to next boundary odqe 
(n(Kl) to B(K2V 


For the remaininq J noints 
outside the boundary find 
the ix)int nearest the edqe 
under consideration? set 
the index number of this 
point eaual to J1 







Fiqure 4b. Block nia<iram of TBIANB, Parts 0 to I. 


Q 


^Ahtnv trum^ 

-ftnu* boumJ.iry odviea?"-^ 
;<*s, sk»n thf ni'xt, 



ConsidcM the 
ary ptunt of 


(I!) 


..id iatv rU boumi- 
the next 

Tf it*s closer to the cui- 
lont ecine than PiJl)# sot 
Jl to its index number (K3) 


r 


T 


( 1 ) 


the 


Consider the adjacent bound- 
ary iK>int or the previous 
edue. It it*s closer to 
adlaccnt edno than P{Jl), 
set Jl to its index number 
(Kl) 



YFS 


-t>l 







Plfiure 4c. Block niacram of TIIANG/ Parts M to o 



(? 

I 

V 

More tho user may insert ad 
(Utionai constraints on the 
triancflo to be formed from 
noint p.Jl. 






• BXock niflorsin 


of TRAINO, Parts R to V 





jO 






Figure 4e. Block Diagram of TRAXN6, Parts X to Y 
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(A) The orocedure beqins with no boundary, no edoes, and all 
points under consideration* Initialize local variables 
and scale the X,Y data. 

J*N, K*L*M»0 
P(j)*j ^or j“l to J 

XMAX»XMIN«XD(1) 

YMAX*YMIN-YD ( 1) 


■"for k*2 to J 
XMAX*maxiroum 
XMIN=minimiun 
YMAXsmaximum 
YMIN-roinimum 


of XMAX,XD(k) 
of XMIN,XD(k) 
of YMAX,YD(k) 
of YMINrYD (k) 


DLXINV=1 . 0/ (XMAX-XMIN) 
DLYINV=1.0/(XMAX-XMIN) 


for k=l to J 
X(k)=(DLXINV) (XD(k) ) 
Y(k)=IDLYINV) (YD(k)) 


(B) 


Begin by taking the last pair of points, (X,Y) in the tisl. 
to be the first boundary point. I'* 


B(1)=J 

J=J-1 

(C) From the remaining points, find the point nearest the first 
il=B(l) 

find i2/il which minimizes 
[;(X(il)-X(i2))" + (Y(il)-Y(i2))‘J 


(D) Now, B(l) to B(i2) is the first edge. There is one edge and 
two boundary points. 

B(2)=i2, K=2, L=l, J=J-1 
E(t,l)=il, E(^,2)=i2 
if i2^J then P(j)=P(j+l) 

for j=i2 to J 
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(E) Now begin circling around the boundary of the polygon^ consid 
ering, in order, each boundary edge. The following Indices 
are maintained — 

K1*B array index of current edge - point 1 
K2-B array index of current edge - point 2 
Bl,B2«index numbers of boundary ooint coordinates 

^ Kl-KT-0 
(y) K1«K1+1 

O K2-K1+1 
^ B1*B(K1) 

B2»B(K2) 

KT»KT+1 

(P) Consider the boundary edge from Bl to B2. For all points not 
yet triangulated (the J points remaining in P) find the point 
that, when triangulated with Bl,B2, minimizes the length of 
the two new edges to be drawn. 


if K1>K then Kl-1 
"n K2>K EHen K2-1 


J1*D1*0 

BFLAG=0 

if J=0 goto (£) 
Tor LJ=1 to J 


PJ=P(LJ) 

then goto 


V' 






^^PJ'^P2^ ‘■^^PJ"^B2^ 


-\ 


® 

(G) 


if Jl=0 or Dl<D then Jl=LJ, Dl.=D 
next LJ 

If less then three edges exist (no triangle defined yet) then 
there are no adjacent boundary points to be considered 


if K;£3 goto (J) 

(H) Consider the adjacent boundary point of the next edge of the 
polygon. Call its index number K3 and see if it's closer 
to the current edge then P(J1). 
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0 j“K3 - K2+1; if K3>K then K3-1 
PK3 » B(K3) 

ii f^PK3"^Bl^ * ^^PK3"^B1^ *^B2*’^Bl!] - 

then goto 

0 ■ -y(XpK3-*Bll‘ * ^ "'pK3'’'b2> 

if J1*0 or D<Dl then Jl»K3» D1»D» BPLAG-1 

closer to the current edge than P(J1) and B(K3) 

0 pK0 = Kl-1; if K0<1 then K0-K 
PK0 = B(K0) 

if r (Vpk0-Ybi) (Xb2-5«Bi) - ‘^PK0-XbiX^B2”^B1) J 1 0-0 


A 


L 


then goto 0 


7 ^ 


D = Y <XpK 0 -^Bl)' ^ (Yp-Yei)' <^B 2 ’^Bi’ ' 

if Jl=0 or D Dl then J1=K0, Dl»D, BFLAG=~1 

(j) Skip the next section if Jl is still zero, since a candidate 
point for trianqulation with edqe Bl,B2 was not found. 

0 if goto {9) 

IKI t£ th« search for “ point^tos 

trEeinr?he^xS'’1ofci;SJe°Shei (J- 0 , . then the „e«t section 
can be ommitted. 

if KT>K or J=0 goto (J) 

^ ^ a.u;eB r^ 4 nt‘ user wfitv insert, an sidditionsl const.r3,int 

(M) At this ^iht the user ^hat one interior an<,le 

be neither very small nor very large. If the 

the test, it is deleted from consideration by setting Jl 0. 

(O) The next procedure checks all boundary edges of the ^lygw 
° fnr intersection with the candidate triangle. If any exist 
ing boundary edge intersects any of the edges to be formed, 
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then the candidate point is rejected. If BFLAG is not zero, 
then the edge defined by J1-K0 or .7l«KJ is exempt from the 
test . 


(N) If there are three or less existing boundary edges or if 
Jl has been set to zero, this test is omitted. 

(ji) if K^3 or J1*0 goto (?) 

if BFLAG»0 then NQ»P(J1) 

IF BFLAG=1 then NQ=B(K3) 

TF BFLAG®-! then NQ=B(K0) 


“for KL=1 to K 
KL=K1 goto 
KN®KL+I; if KL*K then KN®1 



H BFLAGa-l and (KL=K0 or KN®K0) goto 

i£ BFLAG=1 and (KL=K3 or KN=K3) goto 

P1=B(KL) 

P2=B(KN) 


for JL=I to 2 

IT“JL=1 and (BFLAG=0 or BFLAG®!) 
rr JL=2 and (BFLAG®0 or BFLAG®-!) 
if JL=1 then BJ=B! 

FF JL=2 then BJ=2 


XQB=X(NQ)-X(BJ) 

YQB=Y(NQ)-Y(BJ) 

X!2=X(P!)-X(P2) 

Y12=Y(P!)-Y(P2) 

D=XQB »Y!2-YQB*^2 
if 0=0.0 goto 

XIB=X(P1)-X(BJ) 

Y!B=Y(PI)-Y(BJ) 



and KL=K0 goto ^8^ 
and KL=K2 goto ■ 8 





S= (X!B*Y!2-Y!B*X12)/D 
if S<0.0 or S>!.0 goto 




‘^.ext KL 




(?} continue 


(P) If J1 is zero, then the candidate point did not pass the above 
tests or no p^int was found. If BFLAG is not zero, then a 
point on the boundary was found. 


(R) 


if Jl*0 goto ^ 

IT BPLAG^T goto 


_ ■ to ( 

^ BFLAG*- I goto 

The triangulated poTivt is outside the boundary. Establish 
two new edges, a new boundary point and delete one point from 
outside the boundary. 


. 


E(L+1,1) minimum of P(J1), Q(K1) 
E(L+L,2) * maximum of P(J1), B(K1) 
E(L+2,1) ’ minimum of P(J1), B(K2) 
E(L+2,2) = maximum of P(J1), D(K2) 


L=L+2 

KT*0 

M=M+1 

T(M,1) = minimum of 
T(M,2) = middle of 
T(M,3) = maximum of 


P(J1) , B(K1) , B(K2) 
P(J1), B(Kl), B(K2) 
P(J1) , B(K1) , B(K2) 


if Kl^K then B(k+1) = B(k) for k*K to (Kl+1) 


B(K1+1) -=P(J1) 

L«K+1 

J=J-1 


(S) 


il Jl<d~ then P(j)=P(j+l) for j=Jl to J 
goto 


on the boundary . 


trianaulated point is the next point 
EstabUsi; one new edge {from B(Kl) to B(K3)). one new triangle 
(from B(K1) to B(K2) to B(K3)), and delete one point from the 

boundary (B(K2) - 


(^4) KT-^0, KK=0. KKNT-0 

Lwri-i 1) minimum ot B(K1), B(K5) 

( L 1 # 2 ) ma X i mum of B ( K 1 ) ► B l K i ) 
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L*L+1 

K«K-1 

M»M+1 


T(M,1) * minimum of 
T(M,2) “middle of 
T(M,3) * maximum of 


D(K1), B(K2) , B'K3) 
B(K1) . B(K2) . B(K3) 
S(K1) , B(K2) , B<K3) 


if k 2<K then B(k)“B(ktl) for k«K2 to K 
if K2*1 then Kl“Kl+l 


goto @ 


(S) 


The trianqulated point is the previous point on the boundary. 
Establish a new edge (from B(K0) to B(K2)), one new triangle 
Urom B(K0) to B(K1) to u(K21) , and delete a point from the 

boundary (B(K1) . 



KT=0, KK=0, KKNT=0 

E(L+1,1) = minimum of B(K0>, B(K2) 
E(L+1,2) = maximum of B(K0), B(K2) 


L-L+1 

K=KHl 

M=M+1 


T(M,1) = minimum of 

T(M,2) = middle of 
T(M,3) = maximum of 


B(K0) . B(K1) , B(K2) 
B(K0) , B(K1) . B(K2) 
B(K0), B(K1)* B(K2) 


if Kl<K then B(k)=B(k+l) for k=Kl to K 
Kl=Kl-i 


if Kl- L then K1=K 
CP to 


(T) If there are any points remaining outside the boundary, then 

(U) repeat the procedure for the next edge. 

(15) if J>0 and 22^0 ^ 

^ ^ J^O goto 

(V) All points have been triangulated. Check that all boundary 
edges form a concave polygon. 

ijf^ KK^O g oto (55) 
kk«l. kl=o 
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KKNT«KKN'1’+ I 
if KMiT 'j <10tg 




KL-KL+1 

K2»KL<t-l, if K2>K then K2*l 
K1«KL-1, IT Kl- 1 then Kl=K 

PK»B(K), B1»B(K1), B2*B(B21 

ii then 

2£t£ (y) 

if KL<K goto @ 


(X) The tr iiingulet ioti is complete end hes been checked foe e 
concave boundary. Now identify the boundary edges. 



for ^=I to L 
BE ( 5. ) =0 
KL = 0 



KL = KL*-l 

if E{ 1) i‘ B(KL) goto 

KT = Kl+1 

if Kl "K then Kl=l 



H E( • .2) * B(K1) goto 

be; ( .■ ) 
goto 

Kl = KL-l 
i f K 1 * 1 then K I -K 
^ E ( ' , 2 ) f B ( K 1 ) goto 
Mil) = 

if KL' K 




next 



(Y) Finally, establish the indices o' adjacent edges for each 
edge in the triangulation. Each boundary edge will have 
two adjacent edges: each interior, edge will have four. 

initialize TE 

for i = 1 to L 
for i = 1 to 4 

TE ( i) = 0 
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establish TE 
for m»l to M 

for i = to L 

^ E(il/1) ® T(m,l) and Ell,2) * T(m,2) then LI » il 

E(Jl,l) » T(m,2) and E(i,2) » T(ra,3) then L2 * 4 

E(£,l) = T(m,l) and E(t,2) » T(m,3) tEen L3 = t 


X=0; if TE(L1,1)/3 then X»2 

TE(L1,A+1) = L2 
TE(Ll,A+2) = L3 

A*0; ^ TE(L2,1)?<0 then A=2 

TE(L2,A+1) = LI 
TE(L2,A+2) = L3 

A=0; if TE(L3,1)?<0 then A=2 


TE(L3,A+1) = LI 
TE(L3,A+2) = L2 

next m 

RETURN 


Description of Function MIDDLE 


FUNCTION MIDDLE (I,J,K) 

This function is used by the triangulation algorithm to find 
the middle value of throe integer arguments (the value whir-h 
IS neither a minimum or maximum). I,j, and K are assumed 
to be discrete values, no two are equal. 




4.0 


SMOOTHING SUBROUTINE 


The subroutine SMSRF oerforms the optional smoothing of the 

data for the dependent variable. This subroutine uses the 

library routine LLSOF and uses the function POLYX2. 

The smoothing algorithm fits the surface 2 =f(x,y> to a 

poLynoniai c5f the form: 

II K=max(I,J) 

- ^ c-.x^y"' wher.- L=miP (K-i, J) 

i _0 j=o selected parameters 

M 

= :: C;,(xiy3)k 

k=l 

(I+l) (J+l-I/2) J^I 

for M = ' 

(J+1) (I+l-J/2) I>J 

The M terms of the polynomial are each evaluated for n— 1 to N 
noints, where N>M. This evaluation generates an N by M matrix de- 
noted by [AMI. T’he AM matrix is scaled by column so that the mag- 
nitudes of the elements remain close. The scaling factor for each 
column is the average of the absolute values of all elements in 
the column. The N by 1 matrix of z data is known. The task, then, 
is to solve the system 
[AMI [CJ = (Zl 

for the M by 1 matrix C of coefficients. This is accomplished 
by the international Mathematical and statistical Library (IMSL) 
routine LLSQF, which solves the system by means of a linear least- 
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aquarod error criteria. The LLSQF routine is the only Library 
procedure uaod in the contour calculation packaqo. Installations 
which do not have the IMSL library available, would need to replace 
this function with a similar routine. 

After obtaining fCJ from the curve fit subroutine, tho 
coefficients .ire normalized by the same scale factors oriqinally 
used to condition coo*'f icionts are then *-n roolacc 

the orininal z data with new values acquired from evaluation of 
tho polynomial. If the coefficients are not properly found, then 
no smoothing takes place and the oriqinal z data is retained. 

4 . I Description of the Argument List 

CALL SMSRF ( X, Y , Z , ZNKW ,N , I , J , NCOEF,C , IPOIVR, JPOWR) 

Input arguments ; 

X,Y," - arrays containing the function values for z-f(x,y) 

N - tho number of values stored in X,Y,z,:'NEW 

I,J = smoothing parameters selected by the user; used to define 

the K,L values of the polynomial described earlier 

Retu rn Arguments ; 

;'NEW ■ array of new (smoothed) Z data for n=l to N; if the matrix 
" computations fail, ZNEW=Z for all n 

NCOEF = number of coefficients resulting from the values of I and J 

Cj = array of calculated coefficients for i-l to NCOEF 

IPOWRj- for the i=th term of the polynomial, the extjonent of X and 
JPOWR j Y respectively for i=l to NCOEF 



Required Dimensioni: ; 

X(N) 

Y(N) 

Z(N) 

^'^NEW(N) 

B(N) 

AM(C,N) 


IPOWR(C) 

JPOWR(C) 

C(C) 

CNORM(C) 

AVE(C) 

ID(C) 


XX (C) 
H(C) 


4 . 2 D escription of Algorithm 


Fiqures 5a and 5b present a block diagram of the module 
SMSRF. The functions of parts A to J are as follows. 







Figure 


(A) 


(n) 


(C) 


( 0 ) 


(K) 


5a. Block Diagram of SMSRP* Parts A to P 



C START J 


Initialise variables 


^ 

From the I and J esp<'»nont 
values supplied by the user, 
determine the exponent lists 
to be used for the polynomial 


i 

Usimi the exponent lists ami 
the X-Y iiata, construct the 
AM matrix. 




Normali.:e each value in each 
column ot AM by the average 
of that column. 


t 

Set up variables for least 

Siuiartu; 


.. t 

Call tMSl, routine l.LSOF to 
solve (via least squares! 
the system 

|AM|(C1«|2! 
for the matrix C. 


4J 



Piqure Sb. Block Dlaqram of SMSRF« Parts G to J 









(A) Initialize local variables 


<B) 


RN - FLOAT(N) 
if K1 then I«1 
TI J<1 then J-1 
TT « l+l 
J1 » J+1 
NCOEP =0 
lA * 500 

From the I and J exponent values provided by the calling 
program, determine the exponent lists. (IPOWR and JPOWR) 
to be used for the smoothing polynomial. The n-th entry 
in the lists is associated with the n-th term of the 
polynomial. 

K » maximum of II, J1 


I 

L 


for II » 1 to II 
Kll = K-III+1 
L = minimum of Kll, JI 

for JJ = 1 to L 
NCOEF = NCOEF + 1 
IPOWR (NCOEF) = II-l 
JPOWR (NCOEF) = JJ-1 
next JJ 

next II 


(C) Using the exponent lists and the x-y data, construct the 
matrix AM. 


for KCOL = 1 to NCOEF 
lEX = IPOWR (KCOL) 

JEX = JPOWR (KCOL) 


for KROW = 1 to N 
X2 = X(KROW) 

if X2 = 0.0 then X2 = i . 0 
XP a X2**IEX 
Y2 = Y(KROW) 

if Y2 = 0.0 then V2 = 1.0 
YP‘ = Y2**JEX 

AM (KROW, KCOL) » XP*YP 
next KROW 


next KCOL 



Normalize each value in each column of AM by the average 
absolute value of that column. The average of column one 
is always one. 

AVE(l) =1.0 

■for Ll = 2 to NCOEF 

AVE(Ll) = 0.0 

!Tor L 2 * 1 to N , , . , * I 

AVE(Ll) = AVE(Ll) + lAM(L..fLl)l 
AVE(Ll) » AVE(Ll)/RN, if AVE(Ll) » 0, AVE(Ll) 

E or L2 * I to N , 

M(L2,L1) = AM(L2,Ll)/AVE(Ll) 

next Ll 


1.0 


Set up variables for least squares solution. 


M=N 

IER=0 

KBASIS=NCOEF 

tol=o.o 

for KK=1»N B(KK)=?(KK) 

Call IMSL routine LLSQF to 
system AM*C=Z for matrix C 
the solution was found. 


solve (by least squares) the 
If lER is zero on return, then 


CALL LLSQF ( AM , I A , M , NCOEF , B , TOL , KBAS IS , XX , H , I P , I ER) 


if lERj^O then goto 


There were no errors. Transfer 
and divide out the normalization 


the calculated coefficients 
factor. 


Pfor L3 = I, NCOEF 
C(L3) = XX(L3) 

CN0RM(L3) = C(L3) /AVE(L3) 


Evaluate (using function POLVX2) 
a new value for Z. 


each x-y pair to generate 


for L3 = 1 >N 

ZNEW(L3) = -POLYX2 ( 0 ,X(L3) .'t (L3) 


, CNORM , I POWR , JPOWR , NCOEF) 


goto (999 



(J) An error has occurred in the procedure. Set the smoothed 
values of z equal to the original data Set NCOEF to 
negative as an error flag to be checked later. 

ZNBW(i) > Z(l) for t « 1 to N 
NCOEP » -1 
RETURN 


4.3 Description of Function P0LYX2 

FUNCTION P0LYX2 (Z.X.Y.X.IPOWR, JPOWR.N) 

The polynomial evaluation function is used when the smoothing 
option has been invoked. X and Y are the known values of 
the independent variables for which the function value is required. 
Array C is the list of coefficients for each term of the poly» 
nomial. IPOWR and JPOWR are the exponents for each term and N is 
the number of terms. Z is an offset value when evaluating for a 
constant Z. The required dimensions are as follows: 

C(N) 

IPOWR (N) 

JPOWR (N) 



4.4 Description of Subroutine LLSQF 

This is the Library routine taken from IMSL to compute the solution 
of a linear least squares problem. Detailed discussions of the 
argument list and the algorithm can be found in the second volume 
of tne IMSL Library Reference Manual . 
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A summarv of its use is as followsi ^ 

CALL LLSQF (A,IA,M,N,B,T0L,KBASIS,X,H,1P,IER) 


Input Arguments ; 


A 

lA 

M 

N 

B 


roi. 


KBASIS 


M by N coefficient matrix. A is overwritten 
with information generated by LLSQF. 

Row dimension of matrix A as specified in the calling 
program. 

Number of rows in matrices A and B. 

Number of columns in matrix A. 

On input. B is the right hand side of the least squares 
solution (A](X]=[B]. On return, B is overwritten with 
the residual R = B-A*X 

Tolerance parameter to determine the number of columns 
of A to be included in the basis for the least squares 
fit of B. If TOL-0.0 is specified, pivoting is termi- 
nated only if the inclusion of the next column would 
result in a (numerically) rank deficient matrix. 

On input, KBASIS=K implies that the first K columns of 
A are to be forced into the basis. Pivoting is per- 
formea On tne last n-k columns of A. on output:, KBASIS 
gives the number of columns included in the basis. 


Return Arguments ; 


X 

Solution vector 

of length N. 

H 

Work vector of 

length N. 

IP 

Work vector of 

length N. 

lER 

Error parameter 



=0 for normal execution 
= 129 for M<0 or N^^O 
=130 for TOL.-l.O 

(129 and 130 are terminal errors) 
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5.0 


TM>t»gapntAT10N SUBROUTINE 


The subroutine INTERP performs the interpolation of the 
data alona the trlandle edaes. This subroutine use. the 
function POI.VX2 if the smoothinu option has been called. 

The interpolation algorithm is supplied with a set of L ed,e. 
,E(«,l) and EU.2) for t-l to W from the trianpulation. At the 
endpoints of each edge the function value Sj and the independent 
variables and for i-1 to H are known. Additionally, if a 
function has been uenerated tor the values of (from the SMSRF 
subroutine) , the coefficient, and exponents are provided. The 
interpolation procedure will check each edge of the triangulation. 
If the constant value 2 lies between the s function values at the 
end points, then the coordinate. (?j,hj) of 2 relative to the x,y 
coordinates of the endpoints will be calculated, i and n are the 
result of a linear interpolation if the data has not been smooth- 
ed, otherwise, the polynomial previously fitted to the surface is 

solve<^ for the point. 


5.1 


nescrjption of Aroument List 


CALL I NTERP ( X . Y , Z , 2CON s LEDGES , E , I SMOPT , LAMBDA , XI a ETA , J . C , 
IPOWR.JPOWR.NCOEP.M) 
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Input Arguments } 

* the X values of the function Z«f(x»y) 

- the Y values of the function 

- the 2 values of the function 

N » the range of i: the number of points in the X,Y 
and 2 lists 

ZCON ■ the constant value of 2 for which the contour 
values are being interpolated 

LEDGES* the number of triangle edges generated by the 
triangulation procedure 

E(J.»2)* index oointers of endpoints ol each triangle edge» 
e»l to LEDGES 

ISMOPT* smoothing option flag; 1 if SMSRF was called, 0 
if not 

C * coefficients of the polynomial terms as provided 

^ by SMSRF 


Required Dimensions t 


X(N) IE(E,2) 

Y(N) XI (E) 

Z(N) ETA(E) 

LAMBDA (E) 


IPOWR(C) 

JPOWR(C) 

C(C) 


5.2 Description of Algorithm 

Figure 6 presents a block diagram of the module INTERP. 
The functions of parts A to F are as follows: 
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Flc|uir6 6* Block DiacjMin of INTERP 







J-0 


for 2 » 1 to LEDGES 

(A) determine x,y,z values of the endpoints of the next 
edge; order them 


II - E(i,l) , 12 » E(i/2) 


XI 

X2 


X(I1) , 
X(I2) , 


Y1 

Y2 


Y(I1), 
Y(I2) , 


•t 

% 


Z(I1) 

ZU2) 


(B) function values equal or contour value (constant) not 
between endpoints? 



^ Z1 » 22 goto (200 


if Z1 > 22 then reverse XI and X2 

Y1 and Y2 
21 and Z2 

if Z1 > ZCON or ZCON > Z2 goto V2^ 
^ Z2 * ZCON then Z2~= (1 . OOMT) (ZOJfi) 

J = J + 1 

(C) has data been smoothed? 

if not, goto statement label 101 


if_ ISMOPT = 0 goto 


© 


(D) non-linear interpolation is required 
(F) on this edge over the Z surface 

FI = POLYX2 (ZCON, 'Cl,Yl,C,IPOWR,JPOWR,NCOEF) 


for k = 1 to 10 (.1% resolution) 

XN - (Xl+X2)/2.0 
YN = (Yl+Y2)/2.0 

FN » POLYX2( CON,XN,YN,C,IPOWR,JPOWR,NCOEF) 
if FN = 0.0 goto (l^ 

U sign (FI) * sign (FN) then X1=XN, Yl-YN 

if sign (FI) ^ sign (FN) then X2==XN, Y2=YN 

next k 
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- - 1 1 i» Tii-m ■- 



132 

XI (J) - (XH-X2)/2.0 

ETA(J) =* (Yl+Y2)/2.0 



I.AMBDA <• i 

goto 200 


(E) 

(F) 

linear interpolation is required 
on this edge (no smoothing) 


101 

xua, - {-^) - * {^) 

X2 


ETA(J) » ^ Z2-ZT^ * (z2-zl ) 

Y2 


lambda (J) =• ^ 
next n 

return 
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6.0 


CONTOUR SUBROUTINE 


The subroutine CNTOUR draws the reouired contour for z»/S. 

This subroutine calls the user supplied proqram CfJTCRV to draw 
the contour on the qraohics device. 

The contour alqorithm makes use of the results of the 
trianqulation and interpolation procedures in order to establish, 
for each contour to be drawn, the orderinq of the and 
points (for i»l to J) . The coordinates of all intoroolated 
points are known and the trianqulation edqe number associated 
with each coordinate pair is also known. For each edqe, a list of 
adjacent edqe numbers is orovided. A contour line is constructed 
by choosinq a boundary edqe as a startinq ooint (if any) for which 
an interpolated point exists. Then, the remaininq points on the 
contour are ordered by means of searchinq adjacent edqcs for 
interoolated points, until another boundary edqe is encountered. 
For closed contours, the iteration ends if the list of common 
edqes ends. Then a qraphics subroutine is called to draw the 
curve and perform any other user supplied application (for 
example, label the curve). The contour alqorithm then continues 
to the next curve, if there are any points remaininq. This pro- 
cess continues until all contours are drawn and the list of • 
and q coordinates is exhausted. 



6 , 1 Doacription of the Anjuroent List 


CALL CNTOUR ( '/-CON , XI t •i'J’A , LAMBDA , J , I DK , ITE) 
input Arguments ; 


7.C0N 
XI j 

ETAj =“ 

LAMBDA j = 
J 

IBE(t) 
ITK(i,4) ^ 


the constant value of for which the contours 
are beinq drawn 

the x-coordinate of the interpolated point on the 
cdqc C(f)» i LAMBDA! j) 

the y-coordinato of the interpolated ix>int on the 
edge E(e), I =» LAMBDA(j) 

the index number of each edge associated with XI 
and ETA values 

the range of j; the number of interpolated ix>ints 
found for ZCON by the interpolation procedure 

1 if the <.-th edge is a boundary edge; otherwise 
zero 

indices of adjacent edges for the «-th edge 


Kegu i red 0 i im' nj 


ions : 


XI (K) 
KTA(E) 
LAMBDA (F.) 
IBK (E) 
ITE(K,4) 


6 . 2 neseription of Algorithm 


rlgures 7a ami 7e present a block diagram ot the modiile 
CONTOl’R. The functions for parts A to P are as follows. 






Piqure 7b. Block Diaqram of CONTOUR 


Parts II to N 


9 


ntMw tho 4>ftcn contour 
f rom .li to Jr then 
vcixot X 


Are ^ 
there' any 
more tH>intn 
V left; 


9 


y. 

I choose t!u' J-th T'^oint to 


start the cuMitoiir, 



I Seari’h Hu* rom.iinin.i wints? 
— ^ for an a*! iaoont; (common) 

I oiioe. 


-- Was 
an aJjaeH'nt 
eJat' 

V I i>e.i t e\l 


!'ut this i nt eiMH^ a tee! tH>int i 
(M) at the too the list lor , 
this contout ; sot Jl. • 


f — I 

9/ 











(A) Initialize local variable(s) 


J1 * 0 

Search the list of edqcs for a 
boundary edqe. If none is found » 
qo to the procedure for closed contours. 

Jl a J1 I 
LI a LAMBDA (Jl) 

if BE (LI) = 1 qoto (^2' 

^ Jl qoto 

qoto 

(D) Put this point at the top of the list 
and reset Jl . 

© 


Tor JC = l.J 
XT(JC) = xr(jc+i) 

ETA(JC) = ETA(JC+1) 

LAMBDA(JC) = LAMBDA (JCfl) 

JB = J 
L = LI 

(E) Search the remaininq points for an 
adiacont (cominon) odqe. 

JB = JB- 1 

O Jl = 0 

Jl = Jl 1 
Ll - LAMBDAIJI) 

it Ll = TF. (l.,i) (or i = I to 4 qoto © 

^ j 1 J I' goto 

(F) Ah error has occurred. There is no next point. 



^ Jl = J goto (^i) 

XI (J>1) = XI (J) 

ETA(J+1) = ETA(J) 

LAMBDA(J+1)= LAMBDA(J) 



(C) Put this ooint at the top of the list. Continue 
if it*s not a boundary edqe. 




(i) 


XI (J+l) » XI (Jl) 

ETA(J+1) - ETA (.1 if 
LAMBnA(J+l) 3 LAMBDA (Jl) 

for JJ»J1 to J 
XI (.TJ) 3 XI(JJ+1) 

ETA(JJ) 3 ETA(JJ+1) 

LAMBOA(JJ) » LAMBDA(JJ+U 


i£ BE(Ll) y 1 goto 

(H) Draw the open contour from JI 
to J, then reset j. 

NPOIWT = J-JB+1 
ijf NP0INT<1 goto 

Call CNTCRV (XI (JB),ETA(JB) ,NPOINT,7,CON) 



(I) Are there any more points left? 



J = JB-1 
if J<’0 goto 
IT J=0 goto 




(J) Now draw internal linos (closed contours not 
starting or stopping at boundary edges). The 
point at JC = J in the list is chosen to start 
the contour. 




JB = J+1 

(K) Find the next ooint (on the edge with a 
(M) common end point); put it at the top of 

(P) the list; repeat until no more edges are left. 

L = LAMBDA (J) 

JB = JB-l 

Jl = 0/ if JB'J then Jl ^ i 
Jl = Jl^-T" 

LI = LAMBDA (Jl) 

if LI = TE (L.il^r i = 1 to 4, goto ( 14 ) 
IT JKJB qoto VC/ 


(L) Otherwise, no adjacent edge was found; this 
contour.-llne is complete; draw it. 
go to f 17 \ 
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XKJ+l) - XI (Jl) 

ETA(J+l) » ETA(Jl) 
LAMBDA (J+1) a LAMBDA (Jl) 


for JJaJl to J 
XI (JJ) » XKJJ+l) 

ETA(JJ) * ETA(JJ+1) 

LAMBDA (JJ) = LAMBDA (J J+1) 



!■ ^ LI 

If JB ^ 


1 


JJ + JD 

if JB then JJ a JB+l 


(0) Draw the closed contour - the interpolated 
line throuqh the points JJ to J to JJ 


KNT = 0 


to J 


for KK + JJ 
KNT a KNT+1 
XX (KNT) a XI (KK) 
YY(KNT) a ETA(KK) 


NPOINT a KNT+1 
XX(NnoiNT) = XX(l) 
YY (NPOINT) a YYU) 


Call CNTCRV (XX.YY , NPOINT, Z CON) 


(P) Reset J. Establish next contour lines for 
remaininq points or quit if J = 0. 



J a JB-1 
if J''0 goto 

RETURN 









6.3 Description of Subroutine CNTCRV 

This module is supplied by the user and performs the 
graphical nresentation of the contour to the device being used. 
Note that CMTOUR may call this routine several tiroes for each 
constant value of Z^, and a new contour line is provided with 
each call. 


The argument list consists of the following items: 
CALL CNTCRV (XX, YY,NPOINT,ZCON) 

XX = (dimension NPOINT) is the array of X coordinates 
for each point on the contour 

YV = (dimension NPOINT) is the array of v coordinates 

for each point on the contour 

NPOINT = is the number of values provided in the x,y 
coordinate lists 

ZCO'i - is the constant value of Z associated with the 
provided contour line. 
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7.0 PttQGRAMMlWg CONSlDEr »>:10N8 

Th« program* doacribad in thl* doeument have been implemented 
i„ PORTRAN on both an IBM 360/67 (under T8S) and a COC 7600 (under 
SCOPE). The subprogram package, were coded in such a way that a. 
many maehine dependent PORTRAN .tatement. a. possible were elim- 
inated. in lact. the program, appear to be completely portable 
except for (1) the use of IMSL routine bbSCP in SMSRP would need 
to be replaced at installation, where IMSb i. not ..ailable and 
(3) the IBM version uses double precision statements in TBIANG 
that may need modification or deletion. 

The execution time for the contour calculations increases 
with the number of points being processed. The following table 
illustrates typical execution times encountered on a CDC 7600. 

The test cases for this table all made use of the smoothing option 
(With parameters I and J both egual to 2) . and were contrived so 
that three contour lines were generated, each consisting of about 
N/10 interpolated points. The » data points were generated at 
random for these tests. 
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N ® Number of 
data points 

cfx: 7600 
Execution time 
(CPb seconds) 

50 

0.20 

100 

0.47 

150 

0.91 

200 

1.52 

300 

2.66 

400 

4.53 

500 

6.95 


So the execution time is approximately 1.5* (N/200) **1.67 
seconds. 


The algorithms require internal work areas that are used 
to store intermediate calculations durinq execution. The work 
areas required by the triangulation and smoothing procedures are 
the greatest contributing factors to the size of the total object 
time package. The amount of storage required by the trianqulation 
is prop.jrtional to the number of data points to be processed, and 
is approximately equal to 30N. The amount of storaqe required 
by the least-squares curve fitting procedure for smooth data is 
proportional to both the value of N and the maximum number of 
coef ficients to be computed (C) , and is approximately equal to 
C(N+7)+N. The total work area required by all the routines is 
proportional to both C and N, and is approximately N(C+42). 
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For some applications, users may wish to reduce the program 
size. One method, already mentioned, is to eliminate the smooth^ 
inq subroutines if linear interpolation is adequate for the data. 
Size reduction can also be accomplished by decreasing array 
dimensions to accommodate only the maximum number of points and 
coefficients to be processed. Conversely, the array dimensions 
can be enlarged to handle more points and/or coefficients if 
program size is not an imposing consideration. 

Table 1 itemizes all array dimensions which may be given 
new dimensions for the purpose of increasing or decreasing pro- 
gram size as needed. For this table: 

N = Number of datia points to process 
C = Number of coefficients to use 
during smoothing the Z data 
E = 3N-6 = the maximum number of 

triangle edges which can result 
from the triangulation of N points 
T = 2N-5 = the maximum number of 

triangles which can result from 
the triangulation of N points. 


Table 1. Array Dimensions 


Required 

Dimension 


Appears in the 
Following Modules 






Table 2 itemizes local variables that are initialized by 
means of data statements. These data values should be qiven 
new data assignments if any array dimensions are respecified. 


Table 2. Internal Parameter Values 


Data Statement 

Required 

Module Name 

Variable 

Value 


lA 

N 

SMSRF 

MAXCOP 

C 

CNTLNS 

MAXPTS 

N 

j CNTLNS 

1 


As a final note, it should be pointed out that for some 
applications the x and y coordinate values may be used repeatedly 
and only the values of Z will change. For such cases, the x-y 
plane triangulation is valid for each call after the first since 
the triangulation is not based on the Z data. Since the tri- 
angulation can be performed once and then saved, the master 
programs can be easily modified to bypass triangulation of the 
x-y data by inserting an extra parameter in the CNTI.MS argument 
list. Such a scheme would result in a con.siderable savings in 
execution time. 

The subroutine modules described in this report are listed 
in the Appendix. 
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ApnErjnrx 


proc;ram i,rsTiN(\s 




CMTirri3.$» 1 1/05/ 80 09 ; 29x21 


iCjO 

200 


300 

c 

<^00 

c 

500 

c 

600 

c 

700 

c 

BOO 

c 

00 G 

R* 

c 

]000 

c 

HOO 

'L. 

1 200 

c 

1 300 

c 

1^00 

c 

1500 

c 

1600 

c 

imo 

c 

1^00 

c 

1900 

c 

2Q00 

c 

2100 

c 

?700 

c 

2300 

c 

2^00 

c 

2500 

c 

2600 

G 

2700/ 

c 


c 

2900 

c 

3000 

c 

3100 

c 

3 200 

c 

3300 

c 

3^00 

c 

3500 

c 

3600 

3700 

c 


3800 


SUBKuOTlNE lOTLNS (X«Yi Z«Nt 1 SHCPT , IE XP , JEXPt NCNlRSiCLl S I » 

* tPSlOHtIF.RR» 


DRIVER PRCORAM FOR COMPUT IFiG AF:u ORAMInO CCNTGUR LiNEi OF 
CGNSIAtiT L FOP THE FUNCTION I ~ FiXtYI. 


ARv»JHENT LIST OEflNITIUNS - 


X 

Y 

L 

H 

1 SKuPT 
I Exp 
JEXP 
NCNTk:> 

LLIOT 

EP5LUU 

lERK 


= INPUT LIST OF X VALUES 
= INPUT LIST OF Y VALUES 
= INPUT LIST OF I VALUES 

= INPUT SPECIYING THE NUMBER OF VALUES IN A»Y AND x 
= SKUOTHING OPTIuN I LAG CO*NU/OFF, l*YES/uNl 
= I EXPONENT VALUE ICR SHUUTHING 
» J EXPONENT VALUE f OR SMOOTHING 
= NUMBER LF CONTOUR LINES TC BE ORAnN 
IStLF COMPUTING II NCNTRS.LE.Ol 
= LIST 01 CONSTANT CONTOUR VALUES IF NCNTKS.GI.O 

= ERROR Function inlrmalizeu valuei feturneo to 
caller if ismgpt is NoN-ZERO 
= return error flag 

=0 FOR NORMAL RETURN 
* I FOP INVALID VALUE FUR N 

- 2 FOR NUMBER OF ISMGPT COEFFICIENTS **ntATLR 
THAN ‘MAXCOF* OR M 


(NOTE / U NONTPS.LE.O# THEN CLISTIl) * BASE VALJEt 
AND CUST12I » INCREMENT VALUE lOELTAI I 


DlrtENS lo.^ X(N» #Y (NJ »Zlt.) sCL 1ST C 21 
DIMtNSiUN ZNEwCSOOI 





390U . 

AOQO 
4100 
4200 C 
4300 
4400 
4500 C 
4/,00 C 
4 700 C 

4 800 C 
4900 C 

5 GOO C 
5loO 

' ?0C 
5 3C0 
5 400 C 
5 5g 0 C 
5500 C 
5700 C 

5 800 C 
5900 

6 000 c 
61C0 C 
6?oo r 

7.300 c 

6400 

6 500 C 

6t(j0 C 

6 700 f 

6600 C 

69C0 C 

7000 

7100 

7200 

7300 

7400 

75CO 

77,00 


oHL.’ol J 7 

* 

OlMc N 1>1 i./. 


1111494,21.17111494,^1 .All 1494) ,tTA( 14941, 
la M buAl 14941 , I bE 1 14941 

(231 ,JPCJWM23I ,LCEf (231 


l,aTA ha AC of /2 3/ 

data MAAPTS /50C/ 


IMIIAlUE LuCAL VAMAbLES 
ANJ uHfcLA INPUTS FQk EkkORS 


tPSLLN = 0*u 

IF ( f4 • L I • 3 • Ur ♦ T • MAXP T S ) oLTL 


997 


:all iPAiNu lu ifiangjlatf x-y uata p. 

CALL TPlAr*^ I Af Y# NfLtOGtS# It # lt£ t ITL > 

(Cl y 

SMuuTHl’.O kLa^JIHLu? • • 

IF C I iiMuPT, GCTC llO 


( J) 

CHtu*'- ft 


.JtSTLD EXPUNfcNT VALUfcS P CP EKFCrS 


II - lcAr'41 
Jl - ,|rXP4 1 

f«Hl « - M i ( 1 1 » 'I I I 

jUOttll 1 *^c^* c 11AP*11*UEAP41 -IcaP/ 21 
If lil.LKlll OL = (JEAPMl*nLAP.l-JtAP/21 
if iUC.gl.N.'-f •Nt.Gl .MAXCLf 1 oCTO 998 


‘.IS 



C M L Mi »$ I ! 1 / Jb/ 80 09 s 29 : 2 1 


7 700 
7800 
7900 
6000 
a 100 
6 200 

8 300 
6900 
8500 
8600 
8 700 
8800 
3900 
9000 
9100 


C 

C 

C 

C 


OU 12 5 K-1,MAXC0F 
IPDctRlK) = U 
125 JPCim8(K) = 0 

IE I 


1 9200 


Voo 

ZNEkIkI = ZIM 

1 9 300 

G 



1 9900 

C 



1 9500 

C 


(ft 

i 9600 

c 


determine The 

‘ 9 TOO 

c 



96C0 


IdO 

ZMIN ^ zn I 

9900 



ZMAX = ZMIf. 

1 lOOOO 



00 50 K-2,N 

10100 



ZNIN = AMINK 

» !G200 



ZMAX = AMAXll 

] ?030C 


50 

CONTINUE 

10900 

c 



10500 

c 



106 00 

c 


( G • M ) 

10700 

c 


HAS A CLNTUoft 

10300 

c 



10900 



FN » 1 . 0 

11000 



FK = -1,0 

■ moo 


2 00 

IF i NC NTk S , G T , 0 } 

112G0 

c 



1!?00 

c 



1 19U0 

c 


IHI 


call iJ(-RuUTlNE SMSRF TU SMCDTH ThE UATA Z-FCttVI 

CALL 5MSRF < A , Y « Z #i f*EHt M » lEXPi jEXPi NCUEF ,C Lit F i IPOWKt JPiiKHI 
IF INCLEF, I T.OJ CCTO 120 
00 130 K=l,N 

EPiLLM = EPSLON ♦ ( Z ( Ki-ZNEW I K ) ) **2 
1 JO CONTi.iUL 

EPjLON = SORT! EPSLCNI/FLOATIN) 

LGTu 120 

110 00 100 K=l,N 


i>ET ERMINE The RANoE OF THE Z DATA ONOEK CuNSlOcKATl UN 








CU|tNi»lf 11/05/80 09:29:21 

U500 C CALL bJiiKLUTir.E CBVCHK TC VCMF¥ THAT Tht SPtCIHi-J oASfc 

11600 C WALuc IS hITHiN f ANGE Cf DATA, RESET IF NEEOEU 

M 700 C 

11800 call Cb'/CHK (CL I ST 1 1 ) iCCi ST I 21 t ZMlfif ZMAX fCLNEi»l 

11900 IF iCLlSTll ) .Mt.CLNEW) CL I ST ( 1 1 = CLr;E iH 

12000 C 

12100 C (1) 

12200 C JtTEkrtlf.E (NEXTl CUNTUUF CONSTANT VALUE 

’2 300 C 

12A00 210 FK » FK*Uw 

12500 ZCuN « FK*EN*CLIST(2J ♦ CLlST(l) 

12600 IF (ZCUN.OI.ZMlN.Af.U.ZCUN.LT.ZMAXl GOTO 150 

12 700 IF IFN.LT.'J.) GUT G 300 

12800 FK = 0.0 

5 2900 FN = -l.J 

13000 G0T3 210 

13100 C 

rj 13200 l8o X = X*l 

13300 IF iX.GT.NCNlKS) GCTG 300 

13A0G ZCJN * CLISTIK) 

13500 IF (ZCuN.LT .ZMIN.OR.ZCON.GT .ZMAXl UUTU 200 

13600 G 

13700 C Ui 

12800 G CALL SUbROUT INE INTERP TU 

’3900 C INTERPOtATE FOR GGNTUUK LINE GATA POINTS 

lAOOO C 

t^lOO l3J lALL INTEKP IXi Y « ZNEWtNt ZCCNtLEOGESt IE t ISMQPTtLAHBUAt 

14200 • XI rETA.J.CCEF .IPOhR, JPOWRtNCUEFl 

14300 C 

14400 C IK.lI 

’4500 £ ANY DATA POINTS FOUND? • • 

14600 C CALL SUCRUvJT INE CNTQUk TO SuBT THE INTERPGLATEJ PolNTS 

14700 C uM IHE CCNTQUR LINE AND DRAR IT 

14600 C 

14900 SF IJ.Nt.OI GAIL CNTUUR I ZCUN, X I , ET A ,L AMBOA, Jt IBc. ITE I 

1 5000 GOTG 200 

15100 C 

1 5200 ?vja RETURN 




cn TlNi »s f 11/05/ 8 J J9 S 29 : 2 1 


V5300 
15A00 
15500 
1 5600 
15700 


997 I ERR = I 
REIJRN 

998 lERR = 2 
RETJRN 

END 
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SUJKGUTINf: SMSKF ( X » Y i Z i ZNtW tNt 1 1 J i NcQEh tCNuRH. I t J PUhRI 


ICO 

200 C 
300 C 
400 C 
500 C 
400 c 
700 C 
800 C 
900 C 
I CO!) C 
1 100 C 
1200 C 
•300 C 
1400 C 
1500 C 
1600 C 

« *rr,r» ^ 

1 « 

leoo c 
1900 c 
2000 C 

2 i CO C 
2200 C 
2300 C 
2400 C 
2500 C 
2500 C 
2700 C 
28CQ C 
2900 C 
3000 C 
3100 C 
3200 C 
3300 C 
3400 C 
3500 C 
3600 

3 700 
3900 


SUbROUTINL SMSRF PERFCRHS THE GPTIGNAL SMOLTHINC UF UATA BEfUKE 
TklANGULATION OF THE PLANE IS INITIATEO. THE SURFACE DEFINED BY 
Z = FUtY) IS SMOOTHED VIA A POLYNOMIAL CURVE FIT DFFlIUi) BY A 
LEAST SUUAkES CRITERIA. 


ARGJHENTS -- 
UNPUT) 

X,Y,Z ARRAYS OF VALUES DEFINING THE KNOWN SURFACE 

(POINTS IN SPACE FUR THE FUNCTION Z»F(X,YII 

N the number gf points in x,y and Z. 

1,J ARE THE EXPONENTS FOR THE SMOOTHING POLYNuMIAL 

AS SELECTED BY THE USER. 

(RETURN) 

ZNER IS THE ARRAY LF SMOOTHED VALUES FUR THE FUNCTION 

(ZNEW WILL CONTAIN THE ORIGINAL Z UATA UN RETURN 
IF THE SMOOTHING LPERATION FAlLSt IN WHICH CASE 
NCOEF WILL BE SET TO -U. 

NCUEF IS THE NUMBER OF TERMS IN THE POLYNOMIAL RESULTING 

FROM IHE VALUES OF I AND J. NCOEF MUST BE LESS THAN 
Ok EUUAL TO BCTH N AND MAXCOF. 

C IS THE ARRAY OF NCOEF COMPUTED COEFFICIENTS 

IPuwR THE ARRAY OF I EXPONENTS FUR EACH TERM 

JPOWR THE ARRAY OF J EXFCNENTS FOR EACH TERM 

(EACH ELEMENT LF C. IPOWR AND JPuwP IS ASSOCIATED 
WITH THE NCOEF TERMS UF THE POLYNOMIAL# IN uROER) 


DIMEVSIUN X(N) #Y(N) #Z(N) .ZNEWIN) 

DIMENSION )PonR(23) . JPOWR( 2 3 ) # C( 23 ) .CNURMI 23) . AVc( 2 i) 
DlMhNilGN 1P(23) fXX(23) ,H(23) 



SMSPF^t ,n/J5/8iJ 09:29:<|! 


?900 


OIMtN^lGN l5(5oJ) ,AH(5U0|23) 

4000 

c 


4100 


LlAIA lA /5U0/ 

4200 

c 


4300 

c 


4400 

c 

(A) 

4500 

c 

INITIALIZE LLCAL VARIABLES ANU RANGE LHELK 

4600 

c 


4 700 


RfcALN = FLLAT (N) 

4800 


IFll.LT.n I = 1 

4900 


IFIJ.LT.IJ J = 1 

5000 


11 “ l*l 

5100 


Jl ^ J^l 

5200 


NCC6F = J 

5300 

c 


5400 

c 


5500 

c 

IB) 

5600 

c 

DETERMINE THE X AND V EXPONENTS TO BE USED 

5700 

c 

SAVE THEM IN ARRAYS I POhK AND JPQkR 

5800 

c 


5900 


NCUEF = J 

6000 


K = MAXJI II .JI) 

6100 


IF (K.LO.U) GOTO 950 

6200 


DO IdU ll-lill 

6300 


KI 1 = K-I I*l 

6 400 


L NlNOIKIltJl) 

6500 


OIJ 181 JJ-I.L 

6600 


NlllF = NCGEF*'! 

6 700 


IPL'r.P (MLEF) = 11-1 

6800 


JPfJt, K I r^CL EF I * JJ — 1 

6900 


181 CoMlNJE 

7000 


180 CONTINUE 

7100 

c 


7200 

r 


7300 

C 

(C) 

7400 

C 

JSIf.G Thl exponent LISTS FRCiK AbUVE ANJ THE 

7500 

C 

KNOrtN XT DATA PLlNTS, CONSTRUCT Trtt MATRIX AM 

7600 

C 






I 


e 

r 

k 


\ 




7700 



UC 182 l'.LljL=l ,r<LLEf 

7800 



IEa ^ IPuhMKLuL) 

7900 



JtX = jP'jh' 1 KCLL ) 

8000 



OU Za^f KKJW = ltN 

8100 



X2 = XiKKUW) 

0 200 



IF (X2.EU.D.0I X2=1.0 

P300 



XH = X2**1EX 

8^00 



Y2 = YtKKUh) 

esoo 



a lY2.Ey.0.0J Y2=l.O 

6600 



YH = Y2**JEX 

0 700 



AMlKRLiW.KCLL I = XP*YP 

HSOO 


284 

LuMlNUE 

6900 


182 

CUNI INUE 

9000 



KHOrt = l,CCEF 

9100 

c 



9200 

c 



9 300 

c 


(0) 

9A00 

c 


NCiRHALl/E EACH VALUE IN EACH CULUm UF AH BY THF. LLLUHN AVERAGE 

9500 

c 



9600 



AVLIU 1.0 

9 700 



00 t03 LI = 2,MC0EI 

9800 



AVE(Ll) < u.O 

9900 



OU ^02 L2 - 1 iN 

10000 


4U2 

AVEtLlJ = AVE(Ll) ♦ ABS (AHCL2.L1 )J 

lOlOO 



AVE(Ll) ~ AVEILl ) /REALN 

10200 



IF (AVLUn .CU. 0.) AVEILII - l.U 

!0300 



DO A04 L2 = 1 ,N 

10400 


404 

AH(L2iLli AH(L2tLl)/AVElLI) 

1090C 


4U3 

CUNVINUE 

10600 

c 



10 700 

c 



10000 

c 



10900 

c 


(E.F ,G) 

• 1 000 

c 


JSE IHiL ROUTINE LLSUF TO SOLVE IVIA LE AS T-SOUAKE SI 

1 J 100 

c 


THE SYSTEM AM*C = Z FOR MATRIX C 

1 1 200 

c 



moo 

c 



11400 



M = N 
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-j 

-j 


M500 
11600 
1 1700 
1 IRuO 
I ! 900 
12000 
1 2100 
1 2200 
1 2300 
I24C0 
1 2500 
1 2600 
! 2 700 
12800 
12900 
13000 
13100 
13200 
13300 
13A00 
13 500 
13600 
13 700 
13800 
13900 
16000 
16100 
16200 
16300 
16600 
5 6500 
16600 
16700 
I68U0 
1 6900 
150C0 
1 5?00 
15200 


lER = 0 

KBAilS = NCUEF 
TOL = 0.0 


UO 222 KK*1»N 
BIKK) = Z(KK) 

222 CONTINUE 

CALL LLSljF I AMt lAtMiNCOEF ibiTOLt 
IF (ILK.NE.Ol liUTU 950 


KBAS IStXX.H. IP. IE k) 


C 

C 

C 

C 

C 

C 


illlOE LUT THE SCALE FACTOR FROM THE SOLUTION 
uATuiy AMn (^91ABLISH THE COEFFICIENTS 


DU 905 L3 = l.NCOEF 
CIL3) = XXIL3) 

CNURM(L3) = C(L3 )/AVblL31 
905 CONT INUE 


C 

C 

C 

C 

C 

C 


establish THE NEW I VALUES bV 
EVALUATING THE POLYNOMIAL FOR 


EACH KNOWN X-Y PAIR 


DO 936 L3 
ZNEW(L3) 
936 CUNTINUL 
RETURN 


1 tN 

- 1 . 


0*PULYX2( 0.0.X(L3I 


Y(L31 iCNQRM.lPUwRtJPOwKiNCOEFI 


C 

C 

C 

C I J» 

c erkur kltjrn. 

c SENU JACK OLO 

C 

950 00 960 Ll=l.N 
960 ZNEWlLlJ = ZILl) 
UCOLF = -I 


SET NCCEE TO -1 AND 
L VALUES TO CALLING PROGRAM 
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15300 
15400 C 
15500 C 
15600 


RETJRN 

ENO 


CO 
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100 

200 C 
300 C 
^00 C 
900 C 
600 C 
700 C 
ROO C 
900 C 
1000 C 
HOC C 
1200 C 
1300 C 
1900 C 
! 500 C 
1600 C 
1 700 C 
-i 1 800 C 
1900 C 
2000 C 
P’OO C 
2200 C 
2300 C 
2900 C 
2500 C 
2600 C 

2 700 C 
2800 C . 
?90C C 
3000 C 

3 ’00 C 

3 200 C 
3 300 C 
3900 C 
3500 C 
3600 C 
3700 
380D 


SUdKOUTINC TKIANG ( XUt VO » Nt L t E t Bt tT C 1 


9 SET UF N OATA PCINTS ARE KNQHN (X ( I ) t V ( 1 1 1 1* I • U) THEY ARE TO 
BE LONNECTEO BY LINES TU FORM A SET OF TRIANGLES U N.LE. 
MAXPTS). THE FINAL TP I ANGULAT ION ESTABLISHES A CONVEX POLYGON 
OEFlNcD BY LINKED LISTS CF EDGE NUHBERS, ENO POINTS AND 
BOUNDARY EDGES. 


SUBROUTINE INPUT 

XD = ARRAY OF ABSCISSAS 

YO = ARRAY UF ORDINATES 

N = NUMBER UF PUINTS IN X AND Y 

SUBROUTINE OUTPUT 

L = NUMBER OF EDGES LISTED IN Et BE AND TE 
E « LIST OF INDICES OF EACH TRAINGLE EDGE 
BE - 1 IF I OF E IS A BOUNDARY EDGE 

TL > INLICES OF NEIGHBORING EDGES FUR EACH TRAINGlE 
LOCAL VARIABLES 

P = INDICES OF POINTS OUTSIDE THE BOUNDARY 
J NO. OF VALUES IN LIST P 

B - INDEX CF POINTS ON THE BOUNDARY .. INURUEK 
K = NU. UF POINTS LISTED IN ARRAY B 
' T = INOICtS OF ADJACENT TRIANGLE EDGES 
H » NO. UF ROMS USED IN ARRAY T 
X = ARRAY LF SCALED X DATA 

V = ARRAY OF SCALED Y DATA 


IMPLICIT INTEGER (P,BI 
INTEGER T.TE.E 
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2900 C 
AOOO 
A1C0 
4200 
4300 
4400 C 
4 500 C 
4600 

4 700 

4fl00 C 

4900 C 
50U0 C 
5*00 C 
5200 C 

5 300 C 
5400 C 
5500 C 
5600 

c„ 5 700 
o 5800 
5900 
6000 
6100 

6 200 
<300 
6400 
6500 
( 600 
t 700 
6 400 
6^0G 
7000 
7100 
7200 
7300 
7400 
7500 
76Cu 


UlHbNSlu’l 
OIMENSIUN 
OlHENS luU 
DiHLNSION 


a!) 1M • YU (Nil X(500) .Y(500) 
P(5u0) ,b(500) 

E( 1494,21 ,BE( 1494 » ,Tfc ( 1494, 41 
1(995,31 


..UOUBLk PELUSION SPEC IF lOAT ILl. STATizMENTS FCk 1BM36U 
REAC^a T EKM,UC0MP,D,01 ,S,TC 

KEAL*8 XPl ,X?l ,YP1,Y21,XP2,X12,YP2,Y12,XIP,YIP,X2P,Y2P 


(A 1 

THE PKOLEUURL BEoINS WITH NU OCUNOARY, NO EOoES, ANO 

ALL X-Y OAIA POINTS UNDER CcNS lOERA T ION 

SCALl. TriE X.Y UATA ANO irUTIALlZE LOCAL VARIABLES. 


J = N 
K = 0 
L = 0 
M = 0 
KKNT = o 
DO 100 JCNT=l,J 
00 PIJLNTI = JCNT 
XHAX = X j( ! 1 
XMIN - XOdl 
YMAX = YJM I 
YMIN “ YOdl 
JL! 98 K= 2 , N 

XMAX ^ AMAX1(XMAX,XJ(K11 
aMIN = AMIUl (aMIN.XDIKI 1 
YMAX = AMAXII YMAX,YD(Kll 
YMli. = AMINI 1YM(N,Y0{K11 
98 Cu.illNJl 

ULXINV = l.O/(XMAX-XMINl 
UI.YINV = 1, 0/ ( YMAX-YMINl 
JO 96 K-l,f4 
XIK) = aO!K1*DLaINV 
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7700 YCKI = YU(K)*DLYINV 

7600 99 CUNTINUL 

7900 C 

8000 C 

8100 C 

8200 C (8) 

8300 C liEGlN bY TAKING THE LAST PAIR OF POINTS lX.ViJl) IN THE 

8A00 C LIST TO BE THE FIRST BOUNDARY POINT 

8500 C 

8600 B(U = J 

8700 J = J-l 

8800 C 

8900 C 

9000 C 

9100 C (Cl 

9200 C FROM THE REAMINING POINTS. FIND THE POINT NEAREST THE FIRSI 

9300 C 

9A00 C 

oo 9500 12 = 1 

^ 9600 II » Bill 

9700 DMIN = (Xnil-Xilll**2 ♦ I V 1 1 1 l-Y 1 1 II^^Z 

9800 00 270 J1-2.J 

9900 OST = IXC ll I-XIJI n•♦2 ♦ CYI I ll-YCJl ll**2 

’0000 IF ( OST.GE.OMINI GOTO 270 

lOlOO 12 = JI 

10200 UMIN 3 DST 

10300 270 C ONTINUE 

lOAOO C 

10500 C (0) 

10600 C NOW Bin TO BII2I IS THE FIRST EUGE. 

10700 C THEi t IS ONE EDGE AND TtaL BCUNUARY POINTS. 

10800 C 

10900 J » J-l 

llOOC IF II2.GT.JI GOTO 275 

moo 00 276 XNT»I2.J 

11200 P(JCNT) = P(JCNT*1I 

11300 276 CONTINUE 

11600 275 K = 2 
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11500 B(2) = 12 

11600 L = 1 

11700 E(l,l) = HINOIBI 1) .0(2) i 

11600 E(1.2) = MAXOCUd ) ,B(2n 

11900 C 

12000 C 

12100 C 

12200 C (E) 

12300 C N3H BEGIN CIRCLING ARGUNO THE BOUNDARY OF THE POLYGON, 

12A00 C CONSIOEKING, IN ORDER, EACH BOUNDARY EDGE. MAINTAIN THE 

12500 C FULLOulNG INDICES - 

•2600 C K1 = E ARRAY INDEX CF THE CURRENT EDGE - POINT I 

12700 C K2 = B ARRAY INDEX OF THE CURRENT EDGE - POINT 2 

12800 C B1,B2 = INDICES OF BOUNDARY POINT COORDINATES 

12900 C 

13000 Kl' = 0 

13100 XT = 0 

CO 13200 11 Kl = Kl + 1 

13300 IF (Kl.GT.K) Kl*! 

I3A00 12 K2 = Kl*l 

13500 IF (K2.GT.K) K2«l 

13600 Bl = B(K1) 

13700 B2 = BIK2) 

13800 KT - KT+l 

13900 C 

14000 C (F) 

14100 C CONSIDER THE BOUNDARY EDGE FROM Bl TO B2. FDR ALL POINTS NUT 

14200 C TET TRIANGULATED (THE J POINTS REMAINING IN P), FINU THE 

14300 C PCjINl THAT, WHEN TRAINGULATED WITH bl,B2, MINIMIZES T HE LENGTH 

14A00 C jF THE TWO NEW EDGES TO BE DRAWN. 

14 500 r. 

*4600 Ul > 0. 

14700 Jl = 0 

14800 BFLAG = 0 

1490C IF (J.Lw.O) oGTL 6 

15000 00 1 LJ=1,J 

15100 PJ * P(LJ) 

*5200 TERM = ( Y( PJ)-Y ( B I ) I •( X ( B2I -XCU I J)-( X( PJ I-X(BIJ I ♦( Y( B2 ) -Y( Bl 11 
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1 5300 
15400 
1 5 500 
15600 
15 700 
15800 
15900 
160U0 C 
16100 C 
16200 C 
16300 C 
16400 C 
16500 C 
16600 
16700 C 
16800 C 
16900 C 
® 17000 C 

17100 C 
i 7200 C 
17300 C 
1 7900 
17500 
17600 
•77C0 
1 7300 
! 7900 
16000 
1 8100 
18200 
! 3300 
18400 
*8500 C 
18600 C 
18 700 C 
18800 r. 
18900 C 
19000 C 


I 


IF 1 TfcRM.Lt.O. ) GOTO I 

D = SQPTnx(PJ)-x<Bl)M*2*^lYIPj)-yCoin**2) 

IF )*#2) 

IF < Jl .Nt.O. AN0.D1 .lt. D) GOTO 1 

J1 • L J 
01 = D 

continue 


(G) 

IF LESS THAN THREE 
then THERE ARE NJ 
SO Cu TO SECTION J 


EDGES EXIST INO TRIANGLE 
ADJACENT bOUNOARY POINTS 


OfcFlNEO YET), 

TO b£ CONSIDERtO. 


IF (K.LE.3) GOTO 3 


(Hi 


CONSIDER 
POLYGON. 
the CURR 


The ADJACENT BOUNOARY POINT OF THE 
CALL ITS INDEX NUMBER K3 AND SEE 
ENT EDGE THAN P(Jl). 


NEXT EDGE OF THE 
IF ITS CLOSER TO 


K3 « K2+1 
IF 1K3.GT.K) K3=l 
PX3 = 6(K31 

term = IT (PK3 l-Y ( Bl ) )♦ ( X <82 )-Xl bl 

»F (TERM.LE.O.) GOTO 2 
0' SUKTI <X (PK3 )-X<Bl) )**2><Y<PK3 
♦SjRmx<PX3)-X<b2))**2><YIPK3 
I Ji.NE.u.ANO.Ol.LT.O) GOTO 2 
= K3 
» 0 


IF 

J1 

01 


BFLAG = I 


l)-IX<PK3)-X(ttli )*< Ylu2)-YIBin 
)-Y<ail 1 ** 2 | 

I-Ylb2) )♦•2I 


( I 1 

CONSlOEh the 
the POLYGON. 

TO THE CURRENT 


PREVIOUS EDGE OF 
INDEX NUMBER KO ANu SEE IF ITS CLOSER 
EDGE THAN PIJIJ AND B1K3). 



V. i 


I 
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!9!00 
19 200 
’9 300 
19A00 
19500 
19600 
• 9 700 
19800 
19900 
?cnoo 
20100 
20200 
20300 
2CA00 C 
205U0 C 
206C0 C 
20700 C 
20600 C 
20900 
21000 C 
?il00 C 
21200 C 
21300 C 
?U00 C 
21500 C 
2:600 C 
21700 C 
21600 C 
21900 C 
22000 
22100 C 
22200 C 
22300 C 
22AOO C 
2 2 500 C 
22600 C 
2 2 700 C 
22600 


2 CUNT INUE 
Ku = K.l-i 
IF IKu.LT.l ) KU = K 
PKO = ulKO) 

TfcR>t = IV(PK0)-Y(bl)»*(XIB2)-X(bin-(X(HK0J-Aiiiin*( 7(il2»-Y(Bm 
IF I TcRM.Lh .0. ) uCTG 3 

J = S JK r { U(PK0»-XIBU)**2t I Y (PKO)-Y(U 1) l♦♦2) 

2 *SOK ri IX(PK0)-X(B2» I ♦♦2* ( Y ( PKO I -Y (62 ) )**2I 

IF ( Jl .NL.o .A.N(>.Ul.LT.U I CuTU 3 
Jl = KU 
01 = J 
BFt.AG = -1 
j CUI.l INUt 

( J ) 

SKIP iMt rjtXT SECTION If J1 IS STILL IE HL , SINCE A CANDIOATfc 
PCKiT Ftp TklANOULAT ILN WiTh EuGE BltB2 mAS NOT FUUNU. 

IF (Jl.Ej.u) GtTL 9 


(K.L) 

IF THE SEARCH FCH A CANuIOATt PuINT HAS ALkEAUY tUNSluLKLO EACH 
bOUfiOAKV EUGE AT LEAST ONCE (KT.GT.K) Uk IF THE bOJNUARY IS 
uEING CHECKED FOk CONCAVt EDGES IJ=0), THEN THE NEXT SECTlliU 
(SECTILM H) CAN BE OMMITTEO. 

JF ( KT.GT.K. DR. J. EC. U) GOTO 9 

(Ml 

AT THIS POINT THE USER MAY INSERT ANY ADDITIONAL CONSTRAINT 
DN THE TRIANGLE TO BE FORMED BY THE POINT PJl. IF THE 
CANDIDATE TRIANGLE FAILS THE TEST. IT IS DELETED FROM 
CONSIDERATION BY SETTING THE VARIABLE Jl TO ZERO. 

9 CONTINUE 
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I 



22900 C 
23000 C 
23100 C 
23200 C 
23300 C 
23^00 C 
23500 C 
23600 C 
23700 C 
23800 C 
23900 C 
24000 C 
24100 C 
24200 C 
24200 
24400 
24500 
CO 24600 
24 700 
24800 
24900 
25000 
25100 
25200 
25300 
25400 
25500 
25600 
2 5 700 
25800 
25900 
26000 
26100 
26 200 
26300 
26400 
26500 
26600 


(NtUl 

THE NEXT PROCEDURE CHECKS ALL BOUNDARY EDGES OF THE POLYGON 
FOR INTERSECTION HITH THE CANDIDATE TRIANGLE. IF ANY EXISTING 
BOUNDARY EDGE INTERSECTS ANY OF THE EDGES TO BE FORMED BY THE 
CANDIDATE TRIANGLEt THEN THE CANDIDATE POINT IS REJECTED. IF 
BFLAG IS NOT ZERO, THEN THE EDGE DEFINED BY Ji»KO OR Jl=K3 IS 
EXEMPT FORM THIS TEST. 

IF iHtRE ARE THREE OR LESS EXISTING BOUNDARY EDGES UR IF 
Jl HAS BEEN SET TO ZERO, THIS TEST IS OMMITTED. 

IF (K.LE.3.QR.J1.EQ.0) GOTO 7 
IF (BFLAo.EQ.O) NQ <= P(Jl) 

IF (BFLAG.LQ.il NQ » B(K31 
IF (BFLAu.EQ.-I) NU - B(KO) 

DO 108 KCNT-1,K 
IF (KCnT.EU.K1) goto 108 
KS * KCNT+l 
IF (KCNT.EQ.K) KN>1 

IF (BFLAG. EQ.-l. AND. (KCNT.EQ.KO. OR. KN.EQ. KOI) GOTO 108 
IF (BFLAG. EQ. 1 .AND. ( KCNT . EQ.K3.0R .KN. Eu.K3) 1 GOTO 108 
PI * B(KCNT) 

P2 - B(KN) 

DU 8 JLNT«li2 

IF (JCNT.Evi.l .and. (BFLAG. Eg. 0. OR.BFLAG.Eu.il .AND. KCNT. EG. KOI 

• GOTO 108 

IF ( JCNT. EG. 2. AND. I BFLAG. EG. 0. OR. BFLAG. EG. -11 .AN0.KCNT.EU.K21 - 

♦ GOTO 108 
BJ « B1 

IF (JCNT.EG.21 BJ«B2 
XuB = XINQ)-X(BJ) 

YUB = Y(NQ1-Y(BJ) 

X12 » X(Pll-X(P2) 

Y12 = Y(P11-Y(P2I 
0 = XUB*Y12-YQB*X12 




?ftTOO 
?6P00 
?6900 
/70t0 
^'7*00 
?7200 
27?00 
?7^00 
•>7500 
??600 
2 7 700 
?7 6i>0 
?7900 
?8000 
7P!00 
28 200 
? 8?00 
26400 
28500 
28600 
28700 
28800 
28900 
29000 
29100 
29 200 
29300 
29400 
29500 
29600 
29700 
29800 
29900 
30000 
’0100 
30200 
20300 
30400 


C 

( 

C 

C 

C 

C 

C 


n (j.lU.O.) GLTC d 
X Ic - A(Pl)-X(LJ| 

71b = V (Pl)-y (bj) 

j = I Alb*yi2- nB*X12J /I 

17 7'^«bT«0*«LK#S»oT*l»| tiOTC. b 

IC = ( AQn*VlO-YeB*Xlt l/L 

IP i 7C • L T • b * liK • T C # oT • 1 • ) 0( TL b 

Jl - 0 

GLIU 7 

d CJ7.TIf4UL 
»0a LJNTlNJt 
7 LOillNjE 


(Piw) 

»ni/r fUUNO. IF tiFLAC IS NOT ZEku, THEN A 

P01f«r CN The dCUNOAPy mAS fcjuno. 


c 

c 

c 

c 

c 

c 

c 


IF (JI.EJ.U) CLTL 10 
IF (BFlAJ) 150, UO, 4 


THE TRIANoULATEU POINT IS 
Nt« tJOtS, A NEW B0UN0AF.V 
JUTSIJE THE BOUNDARY. 


OUTSIDE THE BUUNOAKy. ESTABLISH I«C 
POINT AND DELETE ONE POINT fkCJM 


*60 EIL*^1,1J 
£IL»1,2) 
EIL«^2 , 1 1 
EI0»2,2) 
KT = J 


HINO(P(Jl),d(Kin 
HAXQ(PUl),B(Kin 
MINOIPIJI ),6CK2) i 
MAXD(P(JU,BIK2J J 


L ^ L*2 
H = Mf-l 

TlM,l) s MI N0( P ( J 1) , B(K1) ,BIK2) ) 
TIM, 21 » MI OULEIP IJI ),bI.M) , bIK2) ) 
TIM, 3) = MAXOIPIJ 1) ,b(Kl) ,b(K2) 1 



I* 1/O‘j/aO J9:i^:0l 


.^0500 If (Nl.L-.r.» LiLTu 190 

■>060: KM ^ r 

ru7no kp: ^ m + i 

•»0fl00 19/ BlrM»l) - LMKHi 

?09U0 K.i - KM-1 

^1000 If (KM»uE»KPl) CitiTi." 197 

•HOT ’ J biM ♦! ) - P ( Jl ) 

’ 1 2uO f. = K ♦ 1 

' ^ ?00 J = J-1 

3 1900 If iJl.ol.J) GUHi lU 

^00 Ou 199 JlNT=JI, J 

iUOO l‘»9 PIJCNT) ^ P(JCNTfl) 

viOTJ lu 

31S00 C 
;1«00 C 

32000 f IS) 

:21C0 r the TKIAMiULATED POINT IS THE NEXT POINT UN THE JOUNOARy. 

222DC C establish ciNE NEh EDGE CFROM B(M) TO BIR3I). uNt NEn 

32300 C TRIANGLE IfKuM BIKl) TO BIK2I TO B(K3))t AND OELcTE ONE POINT 

32900 r fkCH THE UuUNDARY (B(K2)). 

3 2 51.0 C 
3260J C 

•»2700 EILUil) = MlN0IBIK3l#BIKin 

>2600 E(Lfl.2) = HAX0(B(K3)tB(Kin 

22900 K.% = 0 

3 3.)Of' ’ RKNT » 0 

?‘>100 KT = 0 

»3200 L=l*1 

33300 ^=K-l 

? 39 )G M = M* : 

II'Mj - MINJ(DlM)ttJCK2i.B(^3n 
^oOC H M,2) = MI ODLEIb Ir; n ,b(K2) ,EM^3) ) 

■'’/( J U'lfi) = MAa JlblK 1 ) , b(f 2) II (r- :• ) ) 

: I f ( r. ^ . j I .N ) Ot. Tl. 1 5 9 

’3=100 ISl f.CfiT*K2iK 

= 9000 1*51 lU^M) -- ilKLMtl) 

■t.' o (Kj. Ki^Ki-i 

39 300 GUIj 1.J 


T^ I ANL, *$ 


’ l/Jt)/cJ JiC: 


iUVJO 

r 






c 


(RJ 



3 4 500 

c 


!(<;: r-'iANaULAiEu Flint is the pkevilus f'jint un 

THE 

BUUNL'-ARY. 

34?,gO 

L 


:STAol1.H a r.LH EDGE (FRLK BlhOI Tii b(H2)). CNE 

NEm 

TRIANGLE 

34 TOO 

i 


blNu) TL b(Kl» TO oIKZ)), ANi) OELETE uNE 

POINT FROM THE 

34000 

C 


oUo/.jAhY laiKin 



34SOO 

c 





3^»000 


1 5(j 

= HlNOl l(KO) , b(K2) ) 



35100 



EIL»li?) = MAAjlblKO) ,U(K2)) 



35200 



KK - J 



35300 



KKflT = J 



35400 



KT = 0 



35500 



L = L^l 



3 5 600 



K = K-1 



3 5 ?00 



M = rt+1 



35800 



r(M.ll ^ HlNumiKO) >L(Kl) .B(K2)) 



35900 



T(H,2) - M10DLE(b(K0)tB(Kl).b(K2i ) 



360^*0 



riH,3) = HAXOIb(KO) ,6IK1 ) .b(K2) ) 



l(>0 



If IKl.JT.M GlTU 157 



3^ 



UO I5U KCNT»KltK 



3^3uO 


! bd 

BIKCNI) * bIKCta + l) 



36400 


1&7 

K1 » Kl-l 



36500 



IF iKl.LT.il KICK 



?itOO 

c 





H 700 

c 





36000 

c 


i T) 



3e*M0 

r. 


IF THiiil ARE ANY POINTS FEHAlNli.u UUTSIDl THE dOUNUARYt THEN 

3 7000 

c 


KEPLAI IHI PkCCEUURE FOR THE NEXT EDGE. 



3710C 

c 





■»720( 

c 





)7^^j0 


1 j 

IF IJ.lT.G.ANi^.jI.NE.OI GOTO 



^74C J 



IF ij.oi.c) iciu n 



.t 750J 

r 





-7fOi) 

c 





•s 7T-f) 



I J f • 1 rt ) 



37BOC 

c 


ALL POINTS sAVl PLEN T RI ANGOLA 1 LO . CHECK THAI 

Al l 

OUUKOAkV 

»TOC0 

c 


uJv/lS FlRM - ^Li.fAVE POLYOUI*. 




f 







1* : AV, 6»,U/ Jt/do 


■’PlUO 


If- (KK. Ni .01 Ol. Tl 5 5 

rfl’OO 


KK - 1 

’8300 


r 

c 

•P^iJO 

55 

6K9I - 

‘&500 


If* uLTL 17U 

'-.StOO 

5 

K6 = Kl+1 

38 700 


K2 * i^L^l 

:'f<8oo 


ir (K2.uT.<' ) Ki=l 

’8 900 


M - KL-i 

■•9000 


If Ul.U.U Kl-6 

’9100 


PKL » dfKL) 

;9?o0 


bl = 0( K1 ) 

39 300 


is^ B1K2) 

’9900 


T£h4 = «r (fKU-MO: n*lX(B2 J-A161 ) )-(XiPKL l-XIBin*UIb2)- 

39500 


If iTtf y.Ll.O.) CUTI. n 

?9900 


If ( KL.LI .!• ) Guru 5 

39 700 

C 


’9800 

c 


39900 

c 

(XI 

40000 

c 

THE TmI ANUULAT ION IG COMPLETE ANO HAS SEEN CHECKEO FLK 

40100 

r. 

LJNtAVL BOUNOARY, NOb iCLNTIFY THE BOUNDARY EDGES. 

40200 

c 


40 3U0 

C 


40400 

’ 7u 

DO 23 LCNT=l,L 

40500 


bE(LCNT) = 0 

40600 


KL = 0 

'•n 700 

4l 

KL - KL»l 

4080Q 


IF (tlLCNTil J .Nt.BIKLI 1 COTC 22 

4C90C 


» KL*l 

•.1000 


if IKl.CT.KJ Kl«l 

41 iro 


IF (ULCGT *21 .r.t.BCKlM GOTO lt2 

41 200 


BEfLCNlI s 1 

41 =w0 


GOTO 2j 

414u0 

162 

Kl = KL-l 

. : 5.; ,/ 


If IXl.Lr.i) 

•!60G 


If (EILwNI *21 .Gr.BCKin GOTt. 22 

4 1700 


OEUC.M) = 1 

4! 8i;0 


GuT( 2J 








22 

IF (^L.L^.l‘ 1 GOTO 21 



2 3 

..UMIN 

A?1Q0 

(■ 



^2 ?U0 

c 



JOO 

c 


1 Yl 

‘>2'*00 

c 


riNALLt, E5TAPLISH THt Ir.OlUi GF ADJACE.«T EUSES FoK EACH 

h2‘>00 

f 


LOuL lU TMc Tt' lANGULAt IGti. EACH BOUNDARY EDGE hIlL HAVE 


r 


ADJAlLNI edges - EACH INTERIOR EDGE HILL HAVE FUUK. 

4? f 0 

r 



42800 



00 190 Ll =^lt4 

42900 



DU 190 LCNT=1 ,L 

-3000 


i 90 

TEIlCNT iLLI - G 

4 3100 



DC 191 rtwM^ltM 

•• ?2'.'0 



:)U 192 tl^liL 

4 3100 



IF (E(LL •ll.EQ.TIHCNTil).A.NU.EiLLt2)*EQ.TIHCNT«2)l L1>LL 

43400 



IF IE(LL.ll.EQ.T(HCNT»2).A.9D.E(LL.2).EQ.T(MCNT»3n L2=LL 

4 3 40 ' 



IF (EiLL.n.Eu.UHCNT.n. A'4D.E(LLt2).EC.T(HCNT.3n L3»LL 

4 360t 


192 

CONI iNUt 

43 700 



LA'IBDA • 0 

43800 



IF (TElLlil I.NE.O) EAHbDA«? 

43900 



TCILI .LAMBDA^l) > L2 

44000 



TE(Ll(i.AMB0Af2) « L3 

44100 



LAHbOA « 0 

44200 



IF (7EIL2»1 ).NE*01 LAMBDA*2 

44 300 



TE(L2iLAM0UAtl) « Ll 

444G0 



TE(L2tLAMBOA^21 > L3 

44500 



LAMBDA * 0 

♦ 47.00 



IF (TElLitl I.NE.O) LAKBDA - 2 

44700 



TE(L3»LAHBrA»l) = Ll 

44.-,r,o 



TE(l3iLAMBOA»2) » L2 

44900 


1 ->l 

GUN1 INUE 

45 300 

C 



45’Of) 



1 EIJ ’♦ 

45 200 



lnj 



lUOl* »», n/Oi/60 J9:jo:5.'t 


l (»0 
ZOO 
300 
<.C0 
5C0 
MO 
700 
800 
900 
1 000 
’ 100 
1200 
I 300 
! AOO 

I Si!0 

? 700 
1000 
1900 
200C 


fUNCT ICr. MIOOLU I,J,M 
C 
t 

( IfUi. 1-Ji.CTlLN SU[-PkGGPAH IS USEJ BY THfc TRl ANGJLATl 074 ALGORITHH 

lU UNu THt MIDDLE VALUE OF THE TrtlEE INTEGER ARiUMEMS C THE 
( VALUE RrtlCH IS NEITHER A MINIMUM OR A MAXIMUM). It J AND K ARE 

£ ARE ASSJMcU TL til DISCRETE VALUES MlTH NO TWO EUUAL. 

r 

c 


IF (J.LI.I.ANO.I.LT.K) GOTO 100 
IF IK.LT.l.Ar40.1.LT.ji GOTO 100 
IF ( I .LT.J.AND. J.LT.K) GOTO 200 
IF Ik.lT.J.AUj.J.LT.I) GlTU 200 
MIOJLE = K 
kET JRN 

’Oj MIJJLE = 1 
RETURN 

2UU MIUULE 
hETUR.4 
END 


= J 



Zb 


• 1 YX<:»», :i/^t)/<»U J9:3i:l! 


!00 

?00 

5J0 

9C0 

soo 
600 
700 
800 
900 
!000 
1100 
1 ?00 
1300 
190C 
1500 
1600 
1 TOO 
1800 
1900 
?000 
?100 


FUNLTli.N r>LLyX2 ( Z t X » Y • C t 1PCI*8 1 OPCmR»NCCEF I 

f. 

C 

c PCLYX2 IhL POLYNOMIAL EVALUATION FUNCTION USEO WHEN ThE 

C SMQCTH1..0 OPTION HAS BEEN INVOKED. X ANU V LISTS AkE THE 

C KNOhN VALUES OF THE INOEPENUENT VARIABLES. C IS ThE LIST Of 

C COEFFICIENTS FOk EACH TERM. IPOWR AND JPOluR ARE THE EXPONENTS 

C fCi^ EACH TERM ANC N IS THE NUMBER OF TERMS IN THE POLYNOMIAL. 

C Z IS AN OFFSET TERM WHEN EVALUATING FOR A CONSTANT X VALUE. 

C 

C 

UIMENSICN IPGMR (23) (JPONR123) (CC23) 

C 

r 

PULYX2 » J.O 
DO 120 II*ltNCCEF 

P0LYX2 * PLLYX2 ♦ H X«* IPOhK (I I) ) ♦ tV**JPOMRil 1 1 1 1 ♦ CIMl 
120 CONTINUE 

P0LYX2 » Z - PLLYX2 

RETURN 

END 




(, r/C*^^ sSf !l/J5/bv> J9: ?i:17 

IJO SUoKclJ^^^ CbVtHK (ZiLRUt^tLZtilMIfi.ZMAX.ZiNti.) 

2C 0 C 

'vU r 

'-rOO f LUNTuUP cASt VALUt CHtCKlNo FCUTINfc 

500 C **♦*♦*• 

&U0 C 

700 C THIS SUlikOuTINt SHIFTS THE bASE VALUE (ZZERU) UNTIL IT FALLS 

900 C wITmIN The KANuL LiF DATA FCF THIS CUNTOUR II. E. BtTwEEN ZMIN 

900 C ANJ ZMAX). THE SHifTED VALUE ITHE NEW STARTING BASE VALUE I IS 

1000 C RETJPMLD Ti. CALLER AS ZZNEW. THE USER SHIFT INCREMENT COMES 

IVOC C INTJ CbVCMK AS OELZ FOR Z CONTOURS. 

1200 C 

! 300 C arguments - 

1400 C ZZEFO •= oASE VALUE IINPoTI 

I'^OO C OELZ = INCREMENT VALUE (INPUT) 

1600 C ZMIN.iMAX s RANGE OF Z DATA (INPUT) 

1 700 C ZZMh = NEW BASE VALUE* MAY OK ^AY NOT BE 

^ IflOO C THE SAME AS ZZERO (RETURN! 

1900 C 
2000 C 

2100 C — — 

2200 C 
?’O0 C 

2400 IF ( ZMIN. Ew. ZMAX) GOTO 999 

2500 ZZNEH = ZZcRO 

2600 IF (ZHIN.LE.ZZNEH.ANO.ZZNEh.LE.ZMAX) GOTO 999 

2?00 2 IF (ZZNcW.w E.ZHAX) GOTO I 

2RC;. ZZNEtt « ZZNErt ♦ DELZ 

’900 IF (/MIN.IE.ZZNv^. and. ZZNEW. LE. ZMAX) GOTO 999 

3000 GOTC 2 

HOC . 

3iD0 1 ZZNE« = ZZERO 

3300 t IF iZZr.EH.LE. ZMIN! GOTO 999 

•>4U0 ZZNEW * ZZ )Eh - OELZ 

35QG GOTO 999 

: ;3:600 . . :guTC 4 ': 

'■ ri:o c : 

3800 999 RETJKN 







11/05/80 U9:jl:2t 


yQc 

200 

?oo 

C 

9 JO 

C 

5JC 

c 

600 

f 

700 

c 

8G0 

c 

900 

c 

1000 

c 

1100 

c 

1200 

c 

noo 

c 

1400 

c 

’ 500 

c 

1600 

c 

1700 

c 

IBOO 

c 

I960 

c 

2000 

c 

2100 

c 

2200 

c 

2300 

c 

2400 

c 

2500 

c 

2500 

Kr 

2700 

G 

28CJ 

c 

2900 

c 

3000 

c 

. • '0 

c 

3200 

c 

3 •‘CO 

c 

3400 

c 

‘500 

’6C0 

•700 

3300 

c 


SUbRLJTli<f IMfcnP (XtVtUtNtZCONtLEOGEStlEt ISMOPT.LAMOJAtXlt 
» ETAt JfC» IP 3 ^K,JPCiiR,NCaEF) 


SUdftUuTluE iNTERP IS GIVEN A CONSTANT U VALUE (8luUl EOR WHICH 
THE CONTOUR LINE IS TG 8E DRAWN. CHECK ALL GIVEN TRIANGLE EOGESt 
(AKkAY IE) AND CHECK THE VALUES OF U AT THE ENDPOINTS. 

INTERPCLATE FOR ALL POSSIBLE VALUES GN THE TRIANGLE EDGtS. 

IF ISMLPT = 0, then USE A LINEAR INT ERPOLATICNt IF ISMOPT NOT ZERO 
THEN tVALUATE FOR A NON-LINEAR SURFACE USING THE COEFFICIENTS 
FROM SMSRF AND FONCTICN SUBROUTINE POL YX. 

X,Y.u = OEPENLENT AND INOEPEMOENT VALUES FOR 
THE RELATION U»F(X,YI UNPUT) 

4C0N -CONSTANT VALUE OF 2 FOR uHlCH INTERPOLATION 

IS REQUIRED UNPuT) 

LEOGES = MO. CF EDGES IN THE TRIANGOLATION UNPUT) 

IE = EDGE ENDPOINT INDICES FROM TRIANGULATION IINPUTI 

ISMOPT ■ SMCCTHINC OPTION FLAG. 0*CFF* l»ON. UNPUT) 

LAM8DA « INDEX OF EDGES FOR INTERPOLATED POINTS CRETURN) 

XI » LIST OF X-COOROINATES OF INTERPOLATED POINTS 

eta * LIST OF y-COOROINATES OF INTERPOLATED POINTS 

j s NUMBER OF VALUES IN Xlt ETA LISTS 

(Xlt ETA AND j ARE RETURNED) 

C « CIS! OF COEFFICIENTS OF EACH TERM OF THE EQUATION. 

IPUwR.JPUhR are the LIST OF EXPONENTS FOR EACH TERN OF 
THE PCLYNCMIAL JSCD TO SMOOTH THE DATA (INPUT). 

NCULF » NUMfcff OF TERMS IN THE POLYNOMIAL 

tIPCr R,JPCkP,Ct and NCQEF are INPUT) 


DiHEMSiUM X(N) lY IN) lUlM 

UlMENSl- i IE t iA94,2)»XHl-t9HtFTAi 1 9«4) »LANBDAU4VA) 

UiMf.MsiuN i PUkf* { 2>: ) • JPO''>R 1 2 i I » Cl 23 ) 




IMf Rn.$, 


3900 
4 000 
4100 
A 200 
4 300 
A^O 
A5C0 
A6C0 

4 700 
A8G0 
4000 
5000 
5100 
5200 
5300 
5AOO 
55C.G 
5600 

5 ?CU 
f 

5900 
f 300 
fc'OO 
<'200 
eioe 
^400 
f ■?oo 
c400 
700 
rdOO 
^OOO 
7 ’<^0 
7100 
^?00 
730C 
?A00 
75U‘J 
74.00 


e 

c 

c 


^v 

c 

t 

c 


If CN^Ltr.Ll.lJ 

J « 3 


I Si^CPI-O 


wo I I-CNT = ?,l EDGES 

lAi 


OETcRHlNE X,y,2 FOR THE ENOPGINTS OF THE NEXTEDGE - 


oROCR THEN 


11 

12 

XI 

X2 

Tl 

V2 

U1 

U2 


lElLCNT*!) 

lcCLCNT,2) 

XllU 

X112) 

Villi 
Y( 121 
0(11) 

U( 121 


(d) 

FUHCTIU'. 

CONSTANT 


VALc';.S 
IC T 


EwUAL AT ENUPlilNTS CR 
between them? . , 


If 

If 


1 JO 


JEMP 
o2 » 
Ul * 
TENP 
X2 « 
XI » 
TEHP 
V2 » 
VI - 
!f 


IUl.EU.U2) 

1U1.LT.C2) 


GOTC ! 
GCTG 100 


» U2 

Ul 

TENP 
* X2 
XI 

TENP 
» V2 
» VI 

- tenp 

_ .Ul.cF .u^.LT.XCONI GcIu 1 

^ »*02* ' u2 » 1. 000001 * 2C0N 


-) 





' \ 






i- 

>■ 


i 

? * 


] 
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T 7 v.J r 

73tO C 
■^'joo : 
3CJ0 
SIDO C 
8200 f 
8i00 f 
8<f00 C 
8<^00 C 
8600 


rtAS jATA 5EcN SMCCTHfcC? . . 

If NuT, GoTL SECTICN E (STATEHENT LABEL lOll 
If nSHbPI.E^.o) GLIU lOl 


(Difl 

NuN-LlNEAH 1 NT EP PCLAT ICfi IS REgUlREO 
^N T>tIS EUGE UVEF TMt Z-SURfACE 

fl » PlLTXZ »ZCCN»Al»Ylii.tlPCnR»JPU*iR»NCOEFI 


S700 
8 800 
8900 
9000 
9 'CO 
9 200 
9?C0 
9600 
9500 
9600 
9700 
9800 
9900 
lOOCO 
10 ! 00 
’ 0 ?t 0 

• ,300 


00 220 K»l. I J 
XN * |Xl*X2J ♦O.S 
Yf 4 * CYl»Y 21 ♦ 0.5 

FK « PGLYX 2 iZC 0 NtXNt 7 N»C»lPC*<i tJPORR.NCOEF) 
IF tFN.e^.O.I GOTO 132 
IF IFN.lt. 0 ..AND.FI.LT, 0.1 GOTO 235 
IF lfN.GT.u..Afl 0 .Fl. 6 T, 0 .l GOTO 235 
X 2 « XN 
V 2 * YN 
GOTO 220 
235 XI a XN 
VI » YN 
22 v CONTI NO£ 

132 XI Ui » IXl*X2l*0.5 
ETAIJI » IT19Y2I*0.5 
50TG 200 


1 U 6 D 0 C 
Vf 900 C 

IvaOO C U.FI ^ ^ ^ 

’ C ’OC r Ll7«tA»* INTEf PL’LAT ION IS REOUiRtO 

1030C C fGR Trvii EOGi GWEP THE Z-SUPFACi. 



C 

ll /'/ 

u: II » (0..-ALOU) / 1 02-011 


T2 « iZCuN-Ul J/IU2-UII 


•IIJ) a riaAl*T2*X2 

' "»00 

iTAlJJa U*YI+T2*Y;' 

I i 4ii 0 

2JU LA960AI J) = LCNT 








J'. 1 l/yi-mo u9: ^1:?', 


1 !S00 
11600 
: 1 »oo 


Cunt INU5 

fttljRr. 

tNj 
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?00 
^ JO 

c 

r 

SUiiKU^T I/iE 

CwTtuK IZCCN,XI.ETA»LAMBDA, J» IBE* ITEJ 

40C 

V 

c 



§00 

c 

A SET uF J 

INTEKPULATED POK.TS FOR Z-ZCCN (XlllltETAdl ON EDGE 

600 

e 

LAHBDAI 1) 

FOR 1=1«J)» THE CONTOUR LINES MUST NOB BE ORAhN. THERE 

700 

~c 

NAY BE several LlNESi EITHER OPEN uR CLOSED CONTOURS. THIS 

800 

e 

ALGOR I Tn« 

MILL USE THE TRIANGULAT ICN RELAT luNSHlPS TO SORT GUT 

900 

r 

EACH LINE 

IN URDcK. AS LACH CONTOUR LINE IS ESTABLISHED* USER 

1000 

c 

SUPPLIED PROGRAM CNTCRV IS CALLED TO OUTPUT IT TO THE GRAPHICS 

1100 

c 

DEVICE BEING USED. 

1200 

c 



1300 

c 



1^00 

c 



1500 

c 

ARGUMENTS 

(ALL ARE INPUTS) - 

uoo 

c 

ICON 

- Constant value cf i under consi deration 

1700 

c 

XI IJI 

» array of X COORDIANTES OF INTERPOLATED POINTS 

laoo 

c 

ETAU) 

» ARRAY CF Y COORDIANTES OF INTERPOLATED POINTS 

! 900 

c 

LAMBDAlJ) - ARRAY OF EDGE NUMBERS FOR J-TH INTERPOLATED POINT 

2 300 

c 

J 

» HUMBER OF POINTS IN THE LIST OF INTERPOLATED POINTS 

2100 

c 

IBE 

= THE LIST OF BOUNDARY EDGES TAKEN FROM THE TRI ANGULATION 

2200 

G 

JTE 

> LINKED LIST OF INDICES OF ADJACENT EDGES PROVIDED 

2300 

c 


BY THE TRIANGULATIIiN PROCEDURE. 

2400 

c 



2500 

c 

f 



c 

?700 

V 

c 



2«00 

c 



?90Q 


dimension 

X! (I4V4)« ETA (1494) tLAMBuAll 494)* IBEt 14941 ,XX( 14941* 

2000 



Y»(14«4J ,1TE (1494,4) 

3»00 

c 



2 ’00 

c 



3300 

c 



’400 

c 

(AJ 


3500 

r 

IMTIALIZL LOCAL V/RIAslES 

3600 

c 



1700 


IP M.tJ.Oi RtTUR.i 

- 


10 J I = 0 
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o 

o 


3900 f. 
-^000 C 
41QC C 
A?00 C 
4300 
4400 
4500 
4600 
4700 
4900 C 
4900 C 
5000 C 
5100 C 
5?00 C 
5300 C 
5400 
5500 
5600 
5700 
5300 
5900 
6000 
6100 
6200 C 
6300 C 
6400 C 
6500 C 
6600 
6700 
1 300 
6900 
7000 
7100 
7200 
7300 
7400 
7500 C 
7500 £ 


‘eAFCn Ih- Llbl OF tOOES FOF A uCUTiOARY tOGE tbE(il»U 

1 Jl = JUl 

LI = LAMaOAiJll 
IF (l8EtLll.EQ»l> 0170 2 
IF iji.lt.j) goto 1 

sS“cn‘fo« A aoUNOAHV £D« AND PUl 11 A7 THE tOP Of ME LIST. 

PUI MIS INTERPDLAIEO POINT AT THE TOP Of THE 
List FOR THIS CCHTOURi SET Jl 

2 IF (Jl.EU.Ol GCTO 2 
XlUF^U = XlUU 
ETAU>11 » ETAijn 
LAMBDA tJ>U » LAMBDA (Oil 
00 101 JCNT = Jl# J 
XKJCNTI * XltJCNT^U 
ETA(JCNT) » ETACJCNT41I 

lOl LAMBOAtJCNlI = LAMBOAUCNT^l ) 

IUrCH the remaining points for an ADJACENT CCOMMONI EDGE 

3 JIBIG < J 
LCM * Ll 

6 JiaiG * JlBIG-l 
Jl » w 
5 Jl * Jl*l 

kl s LAMaDA(Jl) 

Ot 102 l=l#4 

IF Ui.ci.nElLONT.in GOTO 4 
102 CONTINUE 

NEXT POINT# 
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7700 



IF t Ji.Lr.^^ldlo} GGTC 5 

7800 



6DTU 800 

7900 

C 



3000 

c 


IG) 

8100 

C 


PUT THIS POINT AT THE TOP OF THE 

3200 

c 


LIST. CONTINUE IF ITS NOT A BOUNDARY EDGE. 

8300 

c 



8^00 


4 

xiu^ii = xiun 

8500 



ETAiJ^U « ETAUn 

860G 



LAHB0AU4U » LAMBDA! Jl) 

8700 



DO 103 JCNT » JltJ 

P80Q 



XIUCNT ) * XliJCNT+ll 

8<)0Q 



ETAUCNT) = ETAIJCNT+II 

9000 


103 

LX3dOAUCNT) » LAHbOAUCNT^ll 

«100 



LCNT * LI 

9200 



IF (IBEILU.NE.n GCTL 6 

9300 

G 



^400 

c 



o 9500 

c 


ORAM Tht OPEN CONTOUR LINE THROUGH THE POINTS 

9fc0G 

c 


XiiJUtETAIJU ...... XIIJl+lI fETAtJl*l» ...... 

9700 



THEN RESET J AND CONTINUE 

9800 

c 



9905 

c 



iOOQO 

■ e 



IGIOO 



NPOINT » J-JlBIG+l 

10200 



IF TNPOINT.LE.U GOTO 300 

10300 



CALL GNTORY IXKJIBIGT. ET AC JlBIGIt NPOINT tZCONI 

10400 

c 



10500 

c 



10600 


300 

jaJlBIu - 1 

10700 

( 



10800 

t 


111 

10900 

t 


Afti THEFL ANY MOKE POINTS LEFT? • • 

:!o«)o 

c 



ll!00 



IF t Jl 800t800*10 

11200 

c 



U 3u0 

c 



M4U0 

c 


CJ» 


XlIJItETAUI 




\ .. 


C h T0<> »$ , 1 1 / 0J>/ ao 0 9 : 31 s ft 


11500 
1 1600 
'1700 

iieoo 

11900 

12000 

»2100 

12200 

12300 

12400 

'2500 

12600 

12700 

12800 

12900 

13000 

moo 

13200 
13300 
13400 
12 500 
13600 
13700 
13800 
13900 
140G0 
14100 
14200 
14300 
14400 
14500 


m 
'4 TOC 
148C0 

:49CG 
'!'5 JOf) 
:5ioa 
15200 


C 

c 

c 

c 


NOn JRA». INTLKNAL LINtS ICLOSEO CONTOURS 
jR STJP AT aOUNOAPV EOGESI. THE 

The list li chosen to start 


THE 


That 00 not start 

POINT AT J181G=J IN 
CONTOUR, 


C 

C 

C 

C 

C 

C 


II JlalG = J«l 

LCNI = LAHbOAU) 


PINO THt NEXT POINT PGR THIS CONTOUR ION AN EDGE mITh a rnuMnni 
END POINT). PUT IT AT THE TOP OP THE LIST, AND ReJKt 
NO MORE COMMON EDGES REMAIN POR THIS LINE. 

16 JiaiG » Jlal;,-! 

Jl » 0 

IP IJIBIG.GT.J) Jl^l 
15 Jl . JUl 

U * LAMBDA fJl) 

DO 104 1=1,4 

TP (LI.EO.ITEILCNT,!)) GOTO 14 

104 continue 

IP IJI.LT.JIBIG) GLTG 15 
ID 

OTHERWISE, NO ADJACENT EDGE WAS POUND 

coiJ"!5 ,T. 


C 

C 

c 


14 XIIJ4-1I = XII Jll 
ETAIJ^ll s ETAIJl) 

LAKBOAIJm « LANBOA(Jl) 

DO 105 JCNT » J1,J 
XliJCNT) = XIIJCNT+l) 
ETAIJCNTI * ETALJCNT^l) 

105 LANdOAlJCNTt = i AMbLA f JCNT+l 1 
iCNf « LI 

Tf.TJlB GCTG 16 



11/05/80 09:31:46 
’5300 C 10 J 

15400 G ORAh THl CLOSED CONTOUR LINEi THE INTERPOLATED LINE THROUGH 

’^500 C All Jl»,LTA(\il ) XIUI,ETA(J) ...... XI CJl } •ETAlJll 

15600 C 
15700 C 
15000 C 

15900 C ^ 

16000 i; JJ s JlblG 

l^lOO IE (Jltil3.NE.ll JJ - JlBlG^l 

16200 KNT » 0 

16300 00 510 KK - JJtJ 

16400 KNT * KNT+l 

16500 XXtKNTl « XiiKK) 

16600 VY(KNT) * ETAIKKJ 

16700 510 CONTINUE 

16800 XX(KNT41) - XX(l) 

16900 VY(KNT>11 = VYIU 

17000 NPUINT » KNT*l 

17100 CALL CNTCRV (XXC IltYYd ) .NPUlNTtlCONI 

17200 C — — 

17300 C 
17600 C 

17500 C CP) 

17600 C RESET J. ESTABLISH THE NEXT CONTOUR LINE FOR REMAINING POINTS 

17790 C 3R OiJtT THE PROCEDURE IF NO MORE POINTS REMAIN. 

17800 C 

17900 J s JIBIG - 1 

’8000 IF (J) 800,800,11 

ISIOC 800 RETURN 

5.8200 END 






